Optimizing the output of ode45 by varying a parameter

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

[darova]

Are you asking how to vary R parameter?

[Hayford Azangbebil]

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.

[Hayford Azangbebil]

Massive help, Guru. I really appreciate it.

Answer 1

[Guru Mohanty]

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

[Rachel Ong]

Any idea how I can implmenet this if I have 2 variables to loop?

[Torsten]

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