Deleted vectorized computation feature. Deleted neldermead_contour. Fixed the demos.
[scilab.git] / scilab / modules / optimization / help / en_US / NDcost.xml
1 <?xml version="1.0" encoding="UTF-8"?>
2 <!--
3  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
4  * Copyright (C) 2008 - INRIA
5  * 
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.  The terms
9  * are also available at    
10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  *
12  -->
13 <refentry version="5.0-subset Scilab" xml:id="NDcost" xml:lang="en"
14           xmlns="http://docbook.org/ns/docbook"
15           xmlns:xlink="http://www.w3.org/1999/xlink"
16           xmlns:svg="http://www.w3.org/2000/svg"
17           xmlns:ns4="http://www.w3.org/1999/xhtml"
18           xmlns:mml="http://www.w3.org/1998/Math/MathML"
19           xmlns:db="http://docbook.org/ns/docbook">
20   <info>
21     <pubdate>$LastChangedDate$</pubdate>
22   </info>
23
24   <refnamediv>
25     <refname>NDcost</refname>
26
27     <refpurpose>generic external for optim computing gradient using finite
28     differences</refpurpose>
29   </refnamediv>
30
31   <refsynopsisdiv>
32     <title>Calling Sequence</title>
33
34     <synopsis>[f,g,ind]=NDcost(x,ind,fun,varargin)</synopsis>
35   </refsynopsisdiv>
36
37   <refsection>
38     <title>Parameters</title>
39
40     <variablelist>
41       <varlistentry>
42         <term>x</term>
43
44         <listitem>
45           <para>real vector or matrix</para>
46         </listitem>
47       </varlistentry>
48
49       <varlistentry>
50         <term>ind</term>
51
52         <listitem>
53           <para>integer parameter (see optim)</para>
54         </listitem>
55       </varlistentry>
56
57       <varlistentry>
58         <term>fun</term>
59
60         <listitem>
61           <para>Scilab function with calling sequence
62           <literal>F=fun(x,varargin)</literal> varargin may be use to pass
63           parameters <literal>p1,...pn</literal></para>
64         </listitem>
65       </varlistentry>
66
67       <varlistentry>
68         <term>f</term>
69
70         <listitem>
71           <para>criterion value at point <literal>x</literal> (see
72           optim)</para>
73         </listitem>
74       </varlistentry>
75
76       <varlistentry>
77         <term>g</term>
78
79         <listitem>
80           <para>gradient value at point <literal>x</literal> (see
81           optim)</para>
82         </listitem>
83       </varlistentry>
84     </variablelist>
85   </refsection>
86
87   <refsection>
88     <title>Description</title>
89
90     <para>This function can be used as an external for
91     <literal>optim</literal> to minimize problem where gradient is too
92     complicated to be programmed. only the function <literal>fun</literal>
93     which computes the criterion is required.</para>
94
95     <para>This function should be used as follow:
96     <literal>[f,xopt,gopt]=optim(list(NDcost,fun,p1,...pn),x0,...)</literal></para>
97   </refsection>
98
99   <refsection>
100     <title>Examples</title>
101
102     <programlisting role="example"><![CDATA[ 
103 // example #1 (a simple one)
104 //function to minimize
105 function f=rosenbrock(x,varargin)
106   p=varargin(1)
107   f=1+sum( p*(x(2:$)-x(1:$-1)^2)^2 + (1-x(2:$))^2)
108 endfunction
109
110 x0=[1;2;3;4];
111 [f,xopt,gopt]=optim(list(NDcost,rosenbrock,200),x0)
112
113 // example #2: This example (by Rainer von Seggern) shows a quick (*) way to
114 //             identify the parameters of a linear differential equation with 
115 //             the help of scilab.
116 //             The model is a simple damped (linear) oscillator:
117 //
118 //               x''(t) + c x'(t) + k x(t) = 0 ,
119 // 
120 // and we write it as a system of two differential equations of first
121 // order with y(1) = x, and y(2) = x':
122 //
123 //     dy1/dt = y(2)
124 //     dy2/dt = -c*y(2) -k*y(1).
125 //
126 // We suppose to have m measurements of x (that is y(1)) at different times 
127 // t_obs(1), ..., t_obs(m) called x_obs(1), ..., x_obs(m)  (in this example
128 // these measuresments will be simulated), and we want to identify the parameters
129 // c and k by minimizing the sum of squared errors between x_obs and y1(t_obs,p).
130 // 
131 // (*) This method is not the most efficient but it is easy to implement.
132 // 
133 function dy = DEQ(t,y,p)
134   // The rhs of our first order differential equation system.
135   c =p(1);k=p(2)
136   dy=[y(2);-c*y(2)-k*y(1)]
137 endfunction
138
139 function y=uN(p, t, t0, y0)
140   // Numerical solution obtained with ode. (In this linear case an exact analytic
141   // solution can easily be found, but ode would also work for "any" system.)
142   // Note: the ode output must be an approximation of the solution at
143   //       times given in the vector t=[t(1),...,t($)]  
144   y = ode(y0,t0,t,list(DEQ,p))
145 endfunction
146
147 function r = cost_func(p, t_obs, x_obs, t0, y0) 
148   // This is the function to be minimized, that is the sum of the squared
149   // errors between what gives the model and the measuments.
150   sol = uN(p, t_obs, t0, y0)
151   e = sol(1,:) - x_obs
152   r = sum(e.*e) 
153 endfunction
154
155 // Data
156 y0 = [10;0]; t0 = 0; // Initial conditions y0 for initial time t0. 
157 T = 30;  // Final time for the measurements.
158
159 // Here we simulate experimental data, (from which the parameters
160 // should be identified).
161 pe = [0.2;3];  // Exact parameters
162 m = 80;  t_obs = linspace(t0+2,T,m); // Observation times
163 // Noise: each measurement is supposed to have a (gaussian) random error
164 // of mean 0 and std deviation proportional to the magnitude
165 // of the value (sigma*|x_exact(t_obs(i))|).
166 sigma = 0.1;  
167 y_exact = uN(pe, t_obs, t0, y0);
168 x_obs = y_exact(1,:) + grand(1,m,"nor",0, sigma).*abs(y_exact(1,:));
169
170 // Initial guess parameters
171 p0 = [0.5 ; 5];  
172
173 // The value of the cost function before optimization:
174 cost0 = cost_func(p0, t_obs, x_obs, t0, y0); 
175 mprintf("\n\r The value of the cost function before optimization = %g \n\r",...
176
177 // Solution with optim
178 [costopt,popt]=optim(list(NDcost,cost_func, t_obs, x_obs, t0, y0),p0,...
179                                                        'ar',40,40,1e-3);
180
181 mprintf("\n\r The value of the cost function after optimization  = %g",costopt)
182 mprintf("\n\r The identified values of the parameters: c = %g, k = %g \n\r",...
183                                                                popt(1),popt(2))
184
185 // A small plot:
186 t = linspace(0,T,400);
187 y = uN(popt, t, t0, y0);
188 clf();
189 plot2d(t',y(1,:)',style=5)
190 plot2d(t_obs',x_obs(1,:)',style=-5)
191 legend(["model","measurements"]);
192 xtitle("Least square fit to identify ode parameters")
193  ]]></programlisting>
194   </refsection>
195
196   <refsection>
197     <title>See Also</title>
198
199     <simplelist type="inline">
200       <member><link linkend="optim">optim</link></member>
201
202       <member><link linkend="external">external</link></member>
203
204       <member><link linkend="derivative">derivative</link></member>
205     </simplelist>
206   </refsection>
207 </refentry>