NDcost
======

generic external for optim computing gradient using finite differences



Calling Sequence
~~~~~~~~~~~~~~~~


::

    [f,g,ind]=NDcost(x,ind,fun,varargin)




Arguments
~~~~~~~~~

:x real vector or matrix
: :ind integer parameter (see optim)
: :fun Scilab function with calling sequence `F=fun(x,varargin)`
  varargin may be use to pass parameters `p1,...pn`
: :f criterion value at point `x` (see optim)
: :g gradient value at point `x` (see optim)
:



Description
~~~~~~~~~~~

This function can be used as an external for `optim` to minimize
problem where gradient is too complicated to be programmed. only the
function `fun` which computes the criterion is required.

This function should be used as follow:
`[f,xopt,gopt]=optim(list(NDcost,fun,p1,...pn),x0,...)`



Examples
~~~~~~~~


::

    // example #1 (a simple one)
    //function to minimize
    function f=rosenbrock(x, varargin)
      p=varargin(1)
      f=1+`sum`_( p*(x(2:$)-x(1:$-1)^2)^2 + (1-x(2:$))^2)
    endfunction
    
    x0=[1;2;3;4];
    [f,xopt,gopt]=`optim`_(`list`_(NDcost,rosenbrock,200),x0)
    
    // example #2: This example (by Rainer von Seggern) shows a quick (*) way to
    //             identify the parameters of a linear differential equation with 
    //             the help of scilab.
    //             The model is a simple damped (linear) oscillator:
    //
    //               x''(t) + c x'(t) + k x(t) = 0 ,
    // 
    // and we write it as a system of two differential equations of first
    // order with y(1) = x, and y(2) = x':
    //
    //     dy1/dt = y(2)
    //     dy2/dt = -c*y(2) -k*y(1).
    //
    // We suppose to have m measurements of x (that is y(1)) at different times 
    // t_obs(1), ..., t_obs(m) called x_obs(1), ..., x_obs(m)  (in this example
    // these measuresments will be simulated), and we want to identify the parameters
    // c and k by minimizing the sum of squared errors between x_obs and y1(t_obs,p).
    // 
    // (*) This method is not the most efficient but it is easy to implement.
    // 
    function dy=DEQ(t, y, p)
      // The rhs of our first order differential equation system.
      c =p(1);k=p(2)
      dy=[y(2);-c*y(2)-k*y(1)]
    endfunction
    
    function y=uN(p, t, t0, y0)
      // Numerical solution obtained with ode. (In this linear case an exact analytic
      // solution can easily be found, but ode would also work for "any" system.)
      // Note: the ode output must be an approximation of the solution at
      //       times given in the vector t=[t(1),...,t($)]  
      y = `ode`_(y0,t0,t,`list`_(DEQ,p))
    endfunction
    
    function r=cost_func(p, t_obs, x_obs, t0, y0) 
      // This is the function to be minimized, that is the sum of the squared
      // errors between what gives the model and the measuments.
      sol = uN(p, t_obs, t0, y0)
      e = sol(1,:) - x_obs
      r = `sum`_(e.*e) 
    endfunction
    
    // Data
    y0 = [10;0]; t0 = 0; // Initial conditions y0 for initial time t0. 
    T = 30;  // Final time for the measurements.
    
    // Here we simulate experimental data, (from which the parameters
    // should be identified).
    pe = [0.2;3];  // Exact parameters
    m = 80;  t_obs = `linspace`_(t0+2,T,m); // Observation times
    // Noise: each measurement is supposed to have a (gaussian) random error
    // of mean 0 and std deviation proportional to the magnitude
    // of the value (sigma*|x_exact(t_obs(i))|).
    sigma = 0.1;  
    y_exact = uN(pe, t_obs, t0, y0);
    x_obs = y_exact(1,:) + `grand`_(1,m,"nor",0, sigma).*`abs`_(y_exact(1,:));
    
    // Initial guess parameters
    p0 = [0.5 ; 5];  
    
    // The value of the cost function before optimization:
    cost0 = cost_func(p0, t_obs, x_obs, t0, y0); 
    `mprintf`_("\n\r The value of the cost function before optimization = %g \n\r",...
    
    // Solution with optim
    [costopt,popt]=`optim`_(`list`_(NDcost,cost_func, t_obs, x_obs, t0, y0),p0,...
                                                           'ar',40,40,1e-3);
    
    `mprintf`_("\n\r The value of the cost function after optimization  = %g",costopt)
    `mprintf`_("\n\r The identified values of the parameters: c = %g, k = %g \n\r",...
                                                                   popt(1),popt(2))
    
    // A small plot:
    t = `linspace`_(0,T,400);
    y = uN(popt, t, t0, y0);
    `clf`_();
    `plot2d`_(t',y(1,:)',style=5)
    `plot2d`_(t_obs',x_obs(1,:)',style=-5)
    `legend`_(["model","measurements"]);
    `xtitle`_("Least square fit to identify ode parameters")




See Also
~~~~~~~~


+ `optim`_ non-linear optimization routine
+ `external`_ Scilab Object, external function or routine
+ `derivative`_ approximate derivatives of a function


.. _derivative: derivative.html
.. _optim: optim.html
.. _external: external.html