Using ode45 for a varying input and plot t vs x

Hi,

I am having trouple to plot t vs x with a varying u using the ode45 function. Below is how I attempted to do that. The u given is below

where B= matrix, B^T=transpose of B matrix, Gc(t0,tf)=Gc in the code, x0=initial condition, xf=final condition, psi(tau)=phi_s in the code and psi(tf)=phi_stf in the code

clear all, close all, clc

m = 1; M = 3; l = 2.2; g = 10; d = 0;
t0=0;
tf=2;

syms t
A = [0 1 0 0;0 0 -4.91 0;0 0 0 1;0 0 73.55 0];
B = [0;2;0;-10];
Bt= B';
[M,J] = jordan(A);
phi_s=M * expm(J*t) * inv(M);
Mat = (phi_s*B) * (phi_s*B).';
Gc = int(Mat,0,2)

phi_stf=subs(phi_s, t, tf) % phi_s at final time (tf)
phi_st= % varying phi with time. Having issues making it work

tspan = 0:.001:2;
x0 = [3*pi/4; 0; 0.1; 0]; % initial condition
xf = [0; 0; 0; 0]; % final condition
u=@(x)-Bt*phi_st*inv(Gc)*(x0-phi_stf*xf) % control law

[t,x] = ode45(@(t,x)pend(x,m,M,l,g,u(x)),tspan,x0); % Here is where I am stuck. Phi_st in "u" is supposed to 
% vary with time (tspan) but can't figure out how to add it to it.

plot(t,x(:,1),t,x(:,2),t,x(:,3),t,x(:,4))


function dx = pend(x,m,M,l,g,u)
dx(1,1)=x(2);
dx(2,1)=(u+m*l*((x(4))^2)*sin(x(3))-m*g*sin(x(3))*cos(x(3)))/(M+m*(1-(cos(x(3)))^2));
dx(3,1)=x(4);
dx(4,1)=(-u*cos(x(3))-m*l*((x(4))^2)*sin(x(3))*cos(x(3))+(M+m)*g*sin(x(3)))/(l*(M+m*(1-(cos(x(3)))^2)));
end

Answer 1

[Walter Roberson]

G_c = int(Mat,0,2)

int is always going to return symbolic expression or symbolic function (depending on inputs). But you do not use G_c later in your code

u=@(x)-Bt*phi_st*inv(Gc)*(x0-phi_stf*xf) % control law

That involves the undefined variable Gc (which is a different variable than G_c)

phi_stf=subs(phi_s, t, tf) % phi_s at final time (tf)

The result of subs() is always symbolic expression or symbolic function (depending on the inputs)

You use phi_stf in the definition of u, so if Gc were defined then u would be calculating with a symbolic expression and so would return a symbolic result.

[t,x] = ode45(@(t,x)pend(x,m,M,l,g,u(x)),tspan,x0); % Here is where I am stuck. Phi_st in "u" is supposed to 

and that symbolic result will be passed in to pend as the 6th parameter

function dx = pend(x,m,M,l,g,u)

so inside pend, the u will be symbolic

dx(1,1)=x(2);

The boundary conditions passed by ode45 are always numeric, so x(2) is going to arrive numeric. So dx, which has not been pre-allocated, is going to take the datatype of x(2) which would be double. So dx is going to be double.

dx(2,1)=(u+m*l*((x(4))^2)*sin(x(3))-m*g*sin(x(3))*cos(x(3)))/(M+m*(1-(cos(x(3)))^2));
dx(3,1)=x(4);
dx(4,1)=(-u*cos(x(3))-m*l*((x(4))^2)*sin(x(3))*cos(x(3))+(M+m)*g*sin(x(3)))/(l*(M+m*(1-(cos(x(3)))^2)));

u is symbolic so the right hand side of #2 and #4 are going to be symbolic. Is it going to be possible to convert the symbolic value to double precision? If not then the assignment to dx(2,1) will fail.

It is not good design to knowingly pass in a symbolic number that will result in symbolic calculations taking place that will automatically be converted to double precision because of more or less an accident of the order of assignment to dx. If you want to go ahead with that anyhow, then at the very least pre-allocated dx to force it to double precision, and document what you are doing. But if u is going to be a symbolic number that is certain not to have symbolic variables in it, then why not convert it to double before it reaches pend so that those kinds of issues do not arise?

If Gc happens to be the same as G_c then it will be symbolic and inv(G_c) will be a symbolic matrix inversion... that you are going to perform every time you execute the anonymous function u even though Gc is not changing so inv(Gc) will be constant and should be pre-computed. And you should very likely be using the \ operator instead of multiplying by an inverse matrix anyhow.

u=@(x)-Bt*phi_st*inv(Gc)*(x0-phi_stf*xf) % control law

your phi_st is not defined in your code, so we do not know whether it involves symbolic variables or not.

I suggest that you should be considering using matlabFunction

[mpz]

@Walter Roberson I realized I had a typo for Gc. Should have been Gc not G_c. Also, phi_s is determined above as phi_s=M * expm(J*t) * inv(M); which is symbolic but my issue is how to get phi_st i.e. phi(tau) in the formula of u. Tau is time as it goes from 0 to 2.

So if I understood you correctly, I need to use the matlabFunction to change all my variables that are in symbolic expression before passing u in pend?

[Walter Roberson]

Define phi_st symbolically . Then instead of

u=@(x)-Bt*phi_st*inv(Gc)*(x0-phi_stf*xf)

use

U = -Bt*phi_st*inv(Gc)*(x0-phi_stf*xf)
u = matlabFunction(U, 'vars', {[LIST_OF_VARIABLES]})

where LIST_OF_VARIABLEs might possibly be just t in this situation.