How to optimize the Parameters in my code

Dear All,

I have a problem regarding the optimization of parameters such that the error should be minimized. As given in the code below, I want to find the best values of parameters (UP1 and UP2) in the given range such that my error minimizes. Can anybody help me in this regard. If there is any example similar to my problem, please let me know. I will be very thankful to you.

Best Regards

------------------------------

clear,clc
global XX
tend = 40;   
tspan=1:1:40;

cp1 = 5; cp2 = 3; A1 = 2; rst = 2;  

UP1=.2;  %%%parameter to be optimized range = [.01 1] %%%
UP2=.7;  %%%parameter to be optimized range = [0.5 3] %%%

M =[cp1 0 0 0; 0 cp1 0 0; 0 0 cp2 0; 0 0 0 cp2];

XX=[-rst-cp1*UP1     0       rst         0;
    UP1*cp1          0        0          0;
    rst              0   -rst-UP2*cp2    0;
    0                0      UP2*cp2      0];

 options = odeset('AbsTol',1e-10,'RelTol',1e-10,'InitialStep',1e-2, ...
            'Mass',M,'Jacobian',XX);

u0 = zeros(4,1);
u0(1) = 1;
[t,u] = ode15s(@call, tspan, u0,options);

n1=[u(3,2) u(7,2) u(end,2)];
n2=[u(3,4) u(7,4) u(end,4)];

e1=[0.3  0.8  1];
e2=[0.5  2    5];

%%%% Objective is to minimize the following error %%%%

error = abs((n1(1)-e1(1))^2+(n1(3)-e1(2))^2+(n1(3)-e1(3))^2+...
  (n2(1)-e2(1))^2+(n2(2)-e2(2))^2+(n2(3)-e2(3))^2)

------------------------------

function du = call( t,u ) 
global XX
du=XX*u;
end

------------------------------

Answer 1

[Teja Muppirala]

People often ask how to do this. Basically you need to pass your ODE solver to your optimization solver. This is an example of how to do it:

http://www.mathworks.com/matlabcentral/answers/5054-solving-an-ode-with-best-fit-adjustment-to-empirical-observations

[Teja Muppirala]

If that example doesn't make sense, or you don't know how to apply that to your problem, feel free to ask for more clarification.

[Mathematics]

Dear Teja,

First of all, I am really thankful to you for your kind help. Your example is really good but I could not implement it in my problem. Can you please help me in implementing this in my problem. I have tried the following:

function xout = do_optimization_mine

% The true parameters

e1=[0.3 0.8 1]; % the actual data

T = [3 7 40]'; % need to optimize at three time steps.

hold on;

plot(T,e1);

x0 = [0.2 .7]; % Just some Initial Condition

ub = [5 5]; % Upper bounds

lb = [0 0]; % Lower bounds

F = @(x) COST_mine(x,T,e1);

xout = fmincon(F,x0,[],[],[],[],lb,ub); %<-- FMINCON is the optimizer

legend({'Experimental Data','Fitted Data'});

-------------------------------

function COST = COST_mine(x,T,e1)

A = x(1);

B = x(2);

cp1 = 5; cp2 = 3; A1 = 2; rst = 2;

M =[cp1 0 0 0; 0 cp1 0 0; 0 0 cp2 0; 0 0 0 cp2];

XX=[-rst-cp1*A 0 rst 0;

A*cp1 0 0 0;

rst 0 -rst-B*cp2 0;

0 0 B*cp2 0];

options = odeset('AbsTol',1e-10,'RelTol',1e-10,'InitialStep',1e-2, ...

'Mass',M,'Jacobian',XX);

u0 = zeros(4,1);

u0(1) = 1;

% The cost function calls the ODE solver.

[tout,yout] = ode15s(@call, T, u0,options);

COST = sum((yout - e1).^2);

h = findobj('tag','solution');

set(h,'ydata',yout);

title(['y0 = ' num2str(u0) ' ' 'A = ' num2str(A) ' B = ' num2str(B)]);

drawnow;

-----------------------------------

function du = call( T,XX,A,B )

global XX

du=XX*yout;

end

Answer 2

[Teja Muppirala]

Your answer is actually close, but maybe it wasn't entirely straightforward to match it to the other example without some more experience with seeing these types of problems.

This is (just one) correct implementation of your problem. Make sure you understand why it works. You can see that the answer occurs at a corner, [0.01 3.0].

function [x,f] = dooptim

tspan=1:1:40;

cp1 = 5; cp2 = 3; A1 = 2; rst = 2;  

lb = [0.01 0.5]; %Upper bounds on the variables
ub = [1 3]; %Lower bounds on the variables
x0 = [0.2 0.7] % Some initial condition

mycostfun = @(x) dosolve(x,cp1,cp2,A1,rst,tspan);

[x,f] = fmincon(mycostfun,x0,[],[],[],[],lb,ub,[])

function error = dosolve(x,cp1,cp2,A1,rst,tspan)

UP1=x(1);  %%%parameter to be optimized range = [.01 1] %%%
UP2=x(2);  %%%parameter to be optimized range = [0.5 3] %%%

M =[cp1 0 0 0; 0 cp1 0 0; 0 0 cp2 0; 0 0 0 cp2];
XX=[-rst-cp1*UP1     0       rst         0;
    UP1*cp1          0        0          0;
    rst              0   -rst-UP2*cp2    0;
    0                0      UP2*cp2      0];

options = odeset('AbsTol',1e-10,'RelTol',1e-10,'InitialStep',1e-2, ...
            'Mass',M,'Jacobian',XX);u0 = zeros(4,1);

u0(1) = 1;
call = @(t,u) XX*u;

[t,u] = ode15s(call, tspan, u0,options);
n1=[u(3,2) u(7,2) u(end,2)];
n2=[u(3,4) u(7,4) u(end,4)];
e1=[0.3  0.8  1];
e2=[0.5  2    5];

%%%%Objective is to minimize the following error %%%%
error = abs((n1(1)-e1(1))^2+(n1(3)-e1(2))^2+(n1(3)-e1(3))^2+...
  (n2(1)-e2(1))^2+(n2(2)-e2(2))^2+(n2(3)-e2(3))^2);

[Mathematics]

Dear Teja,

I am very thankful to you for your kind and prompt help.

God bless you.

Best Regards

[yem]

hi i have a system with 2 delays i want to identify these delays with optimisations tools(Newton algorithm)