Optimizing the output of ode45 by varying a parameter

2 次查看（过去 30 天）
Hayford Azangbebil
评论： Torsten 2022-1-7
Dear all,
I'm solving a differential equation using Ode45 and I'm varying a variable R (from 20000 500000) to obtain an optimum (maximum) value of V, an output of the Ode45. I'm, however, at a lost as to how to go about it. Can anyone please offer me any suggestions as to how to go about it?
I will deeply appreciate it.
Thank you.
Hayford.
My code is shown below in the form of a nested function:
function nonlinear_func
Fs=5000;
Ts=1/Fs;
tspan=0:Ts:0.1-Ts; %timespan
R=20000; %load resistance to vary to obtain maximum voltage
IC=[0 0 0]; %% initial condition
%call ode45
[time,state_values]=ode45(@(t,x)nonlinear_func2(t,x,R),tspan,IC);
t=time;
displacement=state_values(:,1); % displacement vector
velocity=state_values(:,2); %velocity vector
V=state_values(:,3); %% generated voltage vector
figure (1)
plot(t,V,'r','linewidth',2);
xlabel('time [s]')
ylabel('Voltage [V]')
title('Generated Voltage Vs Time')
set(gca,'fontsize',12)
figure (2)
plot(time,displacement,'k','linewidth',2);
xlabel('time [s]')
ylabel('Displacement [m]')
title('Displacement Vs Time')
function xdot=nonlinear_func2(t,x,R)
xdot=zeros(1,3);
xdot(1)=x(2);
w=260;
pm=7500;
r=2e-3;
rm=2e-3; %Radius of magnet [m]
hc=6e-3; %height
mt=pi*r^2*hc*pm;
tp=0.25e-3;
lb=38e-3;
b=7.2e-3;
mc=0.0018; %calculated mass
m=mt+mc; %effective masss
Tau1=3415; %yield stress of the fluid
h=0.5e-3; %initial gap between cantilever and magnet in meters
eta=0.288; %viscosity of the fluiding Pa.s
k1=1132; %calculated stiffness of the beam
zeta=0.02; %damping ratio of the cantilever
k2=(4*pi*rm^3*Tau1)./(3*(h+max(x(1)))^2);
c2=(3*pi*eta*rm^4)/(2*(h+min(x(1)))^3);
d31=-315e-12; %d31 coefficient
s11=14.2e-12; %% compliance
phi=0.42;
s11D=s11*(1-phi^2);
z33t=4500; %%relative permittivity constant at constant stress
zeta33=z33t*8.85*10^-12;
z33s=zeta33-(d31^2)/s11; %%permittivity constant at constant strain
Cp=(z33s*lb*b)/tp;
alpha=sqrt((phi^2*Cp*(k1+k2)));
c1=2*zeta*sqrt(m*k1);
wr=sqrt((k1+k2)/m);
Rs=1/(wr*2*Cp);
g=6;
xdot(2)=g*sin(w*t)+(k2*h)/m-(2*alpha*x(3))/m-((c1+c2)*x(2))/m-((k2+k1)*x(1))/m;
xdot(3)=(alpha*x(2)-x(3)/(2*R))*1/(Cp);
xdot=xdot';
end
end
  3 个评论
显示 1更早的评论
Hayford Azangbebil
Hayford Azangbebil 2020-1-9
Yes, Darova. I actually tried with a for loop but the loop wasn't iterating, it actually returns only the last value in the loop. I'm varying R from 20000 to 500000. Within this range, there's a value of R that gives a maximum value of V. I believe using a for loop or any optimization technique will work but I just can't figure out how to do that.
Your help is greatly appreciated.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Guru Mohanty
Guru Mohanty 2020-1-21
Hi, I understand that you are trying to find maximum value of V for a range of R. Here is a sample code. ‘V_max’ gives maximum V and ‘R_val’ gives corresponding R.
clc;
clear all;
Fs=5000;
Ts=1/Fs;
tspan=0:Ts:0.1-Ts; %timespan
R=20000:50000; %load resistance to vary to obtain maximum voltage
IC=[0 0 0]; %% initial condition
V=zeros(length(tspan),length(R));
displacement=zeros(length(tspan),length(R));
velocity=zeros(length(tspan),length(R));
state_values=zeros(length(tspan),3,length(R));
for i=1:length(R)
%call ode45
[time,state_values(:,:,i)]=ode45(@(t,x)nonlinear_func2(t,x,R(i)),tspan,IC);
t=time;
displacement(:,i)=state_values(:,1,i); % displacement vector
velocity(:,i)=state_values(:,2,i); %velocity vector
V(:,i)=state_values(:,3,i); %% generated voltage vector
end
dumm_V=max(V);
[V_max,R_val]=max(dumm_V);
figure (1)
plot(t,V(:,1),'r','linewidth',2);
xlabel('time [s]')
ylabel('Voltage [V]')
title('Generated Voltage Vs Time')
set(gca,'fontsize',12)
figure (2)
plot(time,displacement(:,i),'k','linewidth',2);
xlabel('time [s]')
ylabel('Displacement [m]')
title('Displacement Vs Time')
function xdot=nonlinear_func2(t,x,R)
xdot=zeros(1,3);
xdot(1)=x(2);
w=260;
pm=7500;
r=2e-3;
rm=2e-3; %Radius of magnet [m]
hc=6e-3; %height
mt=pi*r^2*hc*pm;
tp=0.25e-3;
lb=38e-3;
b=7.2e-3;
mc=0.0018; %calculated mass
m=mt+mc; %effective masss
Tau1=3415; %yield stress of the fluid
h=0.5e-3; %initial gap between cantilever and magnet in meters
eta=0.288; %viscosity of the fluiding Pa.s
k1=1132; %calculated stiffness of the beam
zeta=0.02; %damping ratio of the cantilever
k2=(4*pi*rm^3*Tau1)./(3*(h+max(x(1)))^2);
c2=(3*pi*eta*rm^4)/(2*(h+min(x(1)))^3);
d31=-315e-12; %d31 coefficient
s11=14.2e-12; %% compliance
phi=0.42;
s11D=s11*(1-phi^2);
z33t=4500; %%relative permittivity constant at constant stress
zeta33=z33t*8.85*10^-12;
z33s=zeta33-(d31^2)/s11; %%permittivity constant at constant strain
Cp=(z33s*lb*b)/tp;
alpha=sqrt((phi^2*Cp*(k1+k2)));
c1=2*zeta*sqrt(m*k1);
wr=sqrt((k1+k2)/m);
Rs=1/(wr*2*Cp);
g=6;
xdot(2)=g*sin(w*t)+(k2*h)/m-(2*alpha*x(3))/m-((c1+c2)*x(2))/m-((k2+k1)*x(1))/m;
xdot(3)=(alpha*x(2)-x(3)/(2*R))*1/(Cp);
xdot=xdot';
end
  2 个评论
Rachel Ong
Rachel Ong 2022-1-7
Any idea how I can implmenet this if I have 2 variables to loop?
Torsten
Torsten 2022-1-7
A double loop:
for i=1:numel( R)
for j=1:numel(S)
[time,state_values(:,:,i,j)]=ode45(@(t,x)nonlinear_func2(t,x,R(i),S(j)),tspan,IC);
...
end
end

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by