Solving system of ODEs and then get the optimum value for certain parameters.

조회 수: 1 (최근 30일)
prabhat sonu
prabhat sonu 2018년 8월 6일
다시 열림: Walter Roberson 2018년 12월 22일
I want to optimize fro R1 and R2 with the objective function (t_f-1.2e-6).^2+(t_t-50e-6).^2 == 0. Hereis my code that takes R1 and R2 to solve two differential eqns and give me t_f and t_t.
function [t_f,t_t] = fun1(R1,R2)
C1 = 25e-9;
C2 = 1200e-12;
V0 = 10e3;
F = @(t,V) [((1/(R1*C1))*(V(2) - V(1)));((V(1)/(R1*C2))-(V(2)/(R1*C2))-(V(2)/(R2*C2)))];
[t1, V1]=ode45(F,[0 300e-6], [V0 0]);
Vp = max(V1(:,2));
for n=1:size((V1(:,2)),1)
if V1(n,2) == Vp
t_f = t1(n);
elseif ((V1(n,2)-(Vp/2))*1e-3)^2 < 1e-5
t_t =t1(n);
end
end
end
Someone please help me in using optimization tools of matlab like fminsearch or ga to get the optimal values of R1 and R2.

채택된 답변

Eduard Reitmann
Eduard Reitmann 2018년 8월 6일
I am not sure what type of input arguments the "fun1" needs, but it gives the following error when I try fun1(1,1):
Output argument "t_t" (and maybe others) not assigned during call to
"test>fun1".
Both output arguments (t_f & t_t) need to be assigned in the "fun1" function in order to solve the objective function. I made the following changes to "fun1" to assign both outputs (assume zero if not assigned):
if V1(n,2) == Vp
t_f = t1(n);
t_t = 0;
elseif ((V1(n,2)-(Vp/2))*1e-3)^2 < 1e-5
t_t = t1(n);
t_f = 0;
end
I you are happy with the modification above, this code should work theoretically.
fun = @(x) objfun(x(1),x(2));
x0 = [1;1];
options.Display = 'iter';
x = fminsearch(fun,x0,options);
R1 = x(1)
R2 = x(2)
Just remember to define this objective function (after script or separate file, depending on MATLAB version):
function fval = objfun(R1,R2)
[t_f,t_t] = fun1(R1,R2);
fval = (t_f-1.2e-6).^2 + (t_t-50e-6).^2;
end
I got the following answer (you can play with the "options" parameter to increase the accuracy of the answer):
R1 =
1.2505e+03
R2 =
1.3598e+03
  댓글 수: 3
Eduard Reitmann
Eduard Reitmann 2018년 8월 6일
I would rewrite your fun1 code to:
function [t_f,t_t,t,V] = fun1(R1,R2)
C1 = 25e-9;
C2 = 1200e-12;
V0 = 10e3;
F = @(t,V) [((1/(R1*C1))*(V(2) - V(1)));((V(1)/(R1*C2))-(V(2)/(R1*C2))-(V(2)/(R2*C2)))];
[t, V]=ode45(F,[0 300e-6], [V0 0]);
V = V(:,2);
[Vmax,imax] = max(V);
t_f = t(imax);
ihalf = find(V>=Vmax/2,1,'last');
t_t = t(ihalf);
end
The execution code should then be:
fun = @(x) objfun(x(1),x(2));
x0 = [100;100];
options.Display = 'iter';
options.TolX = 1e-20;
x = fminsearch(fun,x0,options);
R1 = x(1)
R2 = x(2)
[t_f,t_t,t,V] = fun1(R1,R2);
figure
plot(t,V,t_f,max(V),'*',t_t,max(V)/2,'*')
t_f
t_t
With the objective function:
function [fval,t,V] = objfun(R1,R2)
[t_f,t_t,t,V] = fun1(R1,R2);
fval = (t_f-1.2e-6).^2 + (t_t-50e-6).^2;
end
Using initial conditions R1=100 and R2=100, I get:
R1 =
2.6172e+03
R2 =
191.7561
This gives a function value of 3.84394e-36 instead of 1.6731e-14 for R1=199 and R2=2.5e3, i.e. more optimized answer.

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Get Started with Optimization Toolbox에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by