Using ODE-45 with control input in state equation updating on each loop

Hello,

I have an issue with my code and I cannot figure out how to solve this issue.

I have a system of equations and the 2 control inputs is constantly updating on each loop. The control inputs have an initial condition of 0.

After the first iteration the control inputs are evaluated using simple algebra of the values of some states.

This is causing me to get all 0's for x(5) and x(6)and these values are needed to calculate u1 and u2 so in return my control values u1 and u2 are also incorrect.

Please help!

function dxdt= odefcn(t,x,u1,u2)

a=4.31e-3;b=1.02e-14;c1=3.41e-10;
f=4.12e-2;g=1.5e-2;h=2.02e1;
k2=6e-1;k3=k2;k1=8e-1;po=2e-11;
s1=1.2e4;s2=7.5e8;delta=1.2e-2;
gamma=9e-1;
r=.5; 



dxdt = zeros(7,1);
dxdt(1)= a*x(1)*(1-(b*x(1)))-(c1*x(4)*x(1))-k3*x(3)*x(1);   
dxdt(2)=-(delta*x(2))-(k2*x(3)*x(2))+s2;
dxdt(3)=-(gamma*x(3))+u2;
dxdt(4)=(g*(x(1)/h+x(1))*x(4))-r*x(4)-(po*x(4)*x(1))-k1*x(4)*x(3)+s1*u1;
dxdt(5)=u1;
dxdt(6)=u2;
dxdt(7)=-1;

%% Run simulatiom
clc
clear
close all

T=60;
tspan = 0:1:T;
x0 = [5e9 1e8 0 1e9 0 0 60]; 

u1=0; u2=0;

 %variable paramaters
 d=.75;duty_cycle=.4;gamma=9e-1;
 
 u1bar=50; u2bar=1;
t=1:1:61;
for k=1:length(t)
     T(k)= (61-k);
     Da(k)=(d*duty_cycle*T(k)*u1bar);
     Db(k)=[d*duty_cycle*T(k)*u2bar];
     D1=transpose(Da);
     D2=transpose(Db);
     u11(k)=(D1(k)-x(k,5)/(gamma*x(k,7)));
     u22(k)=(D2(k)-x(k,6)/(gamma*x(k,7)));
     u1a=transpose(u11);
     u2a=transpose(u22);
     u1(k)= min(u1bar,u1a(k));
     u2(k)=min(u1bar,u2a(k));
end

[t,x]=ode45(@(t,x) odefcn(t,x,u1,u2),tspan,x0);

[Jacob Wood]

It looks like u1 and u2 are both vectors (with length = length(t)), as they were created in the for loop. The expressions for dxdt look like they want u1 and u2 to be scalars instead of vectors. How are u1 and u2 supposed to be calculated?

[T]

Hello Jacob

u1 and u2 should be vectors.

in the code below it shows how they need to be calculated . I think the issue is arising because it uses x5,x6 (look at u11 and u22)

also note that

dxdt(5)=u1;

dxdt(6)=u2;

u1bar=50; u2bar=1;
t=1:1:61;
for k=1:length(t)
     T(k)= (61-k);
     Da(k)=(d*duty_cycle*T(k)*u1bar);
     Db(k)=[d*duty_cycle*T(k)*u2bar];
     D1=transpose(Da);
     D2=transpose(Db);
     u11(k)=(D1(k)-x(k,5)/(gamma*x(k,7)));
     u22(k)=(D2(k)-x(k,6)/(gamma*x(k,7)));
     u1a=transpose(u11);
     u2a=transpose(u22);
     u1(k)= min(u1bar,u1a(k));
     u2(k)=min(u1bar,u2a(k));
end

[darova]

Can you please attach also original equations?

[Jacob Wood]

Mohamed,

I am a little confused on the vector dxdt. You initialize the vector with size = 7x1, but the 5th and 6th elements are both vectors themselves. I think to be compatible with ode45 the dxdt vector needs to be unrolled into a Nx1 vector.

[James Tursa]

@Mohamed: I think you have a misunderstanding of how ode45( ) works. The tspan you pass into it is simply specifying what output points you want returned from ode45( ) ... it does not dictate which points are actually evaluated inside the ode45( ) function. A variety of calls to your derivative function at times other than tspan will be called. So your derivative function needs to be able to calculate or use controls for an arbitrary time that ode45( ) will pass into it, and not just at the tspan points. You need to modify you derivative function to operate in this manner. Maybe that means interpolating your u vectors for a specific time, or calculating u values on the fly based on the time, etc. You will need to decide this and code it appropriately.

Answer 1

[T]

I figured it out with the help of a collegue

I calculate the control within the function code

also I treat the controls as consts

D1=d*duty_cycle*T*u1bar;
u11=(D1-x(5)/(gamma*x(7)));
u1max= min(u1bar,u11);

[darova]

They also should be inside of ODE function

[darova]

[T]

thank you Darova!

Answer 2

[T]

Attached are the equations

[darova]

Can you explain more?

What is the difference between , and ?

[T]

assume u1= u1max and u2=u2max and u1bar,u2bar are constans