Integration of a vector inside function for ode45

I am trying to solve a set of differential equations using ode45 command. The problem is I want to integrate one of the variables of the differential equations. For that I need to save the values/solution of that variable in a vector inside the function. I am unable to store x(1) as a vector. In the code, I need to integrate x(1) for which I have used cumtrapz but its not helping . My differential equations arre x1dot=x2 ;x2dot=-9.81*sin(x1)+u where dot means diffeentiation with respect to time. For u, e, ei and de are as defined in the code. Kindly help.

function dxdt=slmc(t,x)
dx1dt=x(2);
e=-x(1);
de=-x(2);
ei=0.01*cumtrapz(e);
u=75*e+5*ei+25*de;
dx2dt=-9.81*sin(x(1,:))+u;
dxdt=[dx1dt;dx2dt];
end


t=0:0.01:10;
in=pi/2;
indxdt=0;
[t,x]= ode45(@(t,x) slmc(t,x),t,[in indxdt]);
plot(t,x(:,1))

Answer 1

[Torsten]

In the code below, x(3) = integral_{tau=0}^{tau=t} -0.01*x(1) dtau

t=0:0.01:10;
in=pi/2;
indxdt=0;
[t,x]= ode45(@(t,x) slmc(t,x),t,[in indxdt 0]);
plot(t,x(:,3))

function dxdt=slmc(t,x)
dx1dt=x(2);
e=-x(1);
de=-x(2);
ei = x(3);
u=75*e+5*ei+25*de;
dx2dt=-9.81*sin(x(1,:))+u;
dxdt=[dx1dt;dx2dt;-0.01*x(1)];
end

[Paul]

Should dxdt(3) be +0.01*x(1) ?

[Torsten]

I think -0.01*x(1) should be integrated to get ei due to the OP's setting

ei=0.01*cumtrapz(-x(1));

, shouldn't it ?

[Paul]

Correct. My mistake.

[SM]

I want to integrate e=-x(1) with a time step of 0.01. Could you clarify why another variable x(3) needs to be introduced for this? This e should be a vector otherwise the integration yields to be zero at each instant.

[SM]

I want to see the response of x(1). So why in the plot command x(:,3) is witten?

Thank you for your inputs.

[Torsten]

Because out of interest, I chose x(3) to be plotted.

Be uninhibited - it's your code. Just change x(:,3) to x(:,1).

[Torsten]

I want to integrate e=-x(1) with a time step of 0.01. Could you clarify why another variable x(3) needs to be introduced for this?

You want to integrate e = -x(1). Thus you want to get

E(t) = integral_{tau=0}^{tau=t} e(tau) dtau

If you differentiate E(t) with respect to t, you get (by the fundamental theorem of calculus)

dE/dt = e(t) = -x(1).

So you have to solve a third differential equation

dx(3)/dt = -x(1)

where E is set to x(3).

Since in your code, you chose ei = 0.01*cumtrapz(-x(1)), I integrated -0.01*x(1) instead of -x(1).

But I'm sure that after my explanation, you will be able to change this easily on your own (although it's not a mistake).

[SM]

Thanks a lot @Torsten , I undestood now. However I still have some doubts. Like I want to verify whether what I intended to do with my code, is correct or not. That 'u' specifies a PID controller for the second order system as given by the differential equations mentioned in the question and the code. When I try to obtain results in an alternative way, without using ode45, I am getting different results. Could you help me with this?

I have taken an example problem for this:

Without using ode45:
A=[-10 -20;1 0];
B=[1;0];
C=[0 1];
D=0;
Kp=50;Ki=30;Kd=16;
c1=pid(Kp,Ki,Kd);
G=ss(A,B,C,D);
% initial(G,[pi/2 0])
closed=feedback(G*c1,1);
initial(G,[pi/2 0])

With using ode45:
function [dxdt,y]=test_pid2(t,x)
e=-x(1)
de=-x(2);
ei=x(3);
u=50*e+30*ei+16*de;
dx1dt=-10*x(1)-20*x(2)+u;
dx2dt=x(1);
dxdt=[dx1dt;dx2dt;-x(1)];
y=[0 1]*[x(1);x(2)];
end

t=0:0.01:10;
in=pi/2;
indxdt=0;
[t,x]= ode45(@(t,x) test_pid2(t,x),t,[in indxdt 0]);
[~,y]=cellfun(@(t,x) test_pid2(t,x.'),num2cell(t),num2cell(x,2),'uni',0);
r=zeros(size(t));
yi=cell2mat(y)
figure(2)
plot(t,x(:,1),'b')
figure(3)
plot(t,x(:,2),'k')
figure(4)
plot(t,yi,'g')

[Torsten]

I don't have the time to look for what pid, ss and feedback do.

But I can list the ODE system you are trying to solve. Maybe you can find the mistake therein.

You solve

dx1/dt = -10*x1 - 20*x2 - 50*x1 - 30*(integral_{tau=0}^{tau=t} x1(tau) dtau) - 16*x2
dx2/dt = x1

with initial conditions

x1(0) = pi/2
x2(0) = 0

My guess is that it is incorrect.

[Paul]

The initial() command is plotting the IC response of the plant, not the closed loop system.

The closed loop system formed by G and c1 is not the same as that implemented in test_pid2. In the latter, the control input is

u=50*e+30*ei+16*de = -(50*x1 + 30*intx1 + 16*x2) % that 16*x2 was probably meant to be 16*x1dot

whereas in the former the control is

u = -(50*y + 30*inty + 16*ydot) = -(50*x2 + 30*intx2 + 16*x2dot)

[SM]

Thanks @Torsten and @Paul for your valuable inputs!