How to get the state trajectory with a GIVEN input VECTOR using ode45 solver
13 views (last 30 days)
I would like to check whether my previous trajectory is plausible or not.
From the previous calculations I obtained a 50x1 double VECTOR I called u_opt , which represents the input (u) of the system of equations You see below.
function f = calcf(t,x,u,param)
% Physical parameter in the struct "param"
M = param.M;
g = param.g;
L = param.L;
m = param.m;
d = param.d;
Sx = sin(x(3)); % sin(x)
Cx = cos(x(3)); % cos(x)
N = M + m*(1-Cx^2); % denominator
f = [x(2); (1/N)*(-m*g*Cx*Sx + m*L*x(4)^2*Sx - d*x(2)+ u); x(4); (1/(N*L))*((m+M)*g*Sx -Cx*(m*L*x(4)^2*Sx - d*x(2)) - Cx*u)];
Now I would like to hand this vector u_opt over to the ode45 solver to obtain a 50x4 double x_ode that I can compare with my previous trajectory and see if it's plausible or not as follows:
[t_ode,x_ode] = ode45(@(t,x)calcf(t,x,u_opt,param),tspan,x0);
with tspan and x0 being just
tspan = linspace(0,tf,N); % N=50 - tf=5s ... final time
x0 = [0;0;pi;0]; % start state - x0_1 position , x0_2 velocity , x0_3 angle , x0_4 angular velocity
The Command window shows this error
Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error using odearguments (line 95)
@(T,X)CALCF(T,X,U_OPT,PARAM)returns a vector of length 102, but the length of initial conditions vector is 4. The vector returned by
@(T,X)CALCF(T,X,U_OPT,PARAM) and the initial conditions vector must have the same number of elements.
The ode45 solver seems to work with a given input being a scalar, but I would like to calculate the solution of the diff. equation for every state at a time 't' being element of tspan, so that I ca compare my results / plots. That's why I want the vector u_opt to be the input.
Can anyone help me?
Thanks and cheers!
Ayush Gupta on 10 Sep 2020
ode45() is not suitable for finding multiple solutions, except through the mechanism of running multiple times with different u until the number of distinct solutions found has accumulated to the number desired. A loop can be run over the ode and the solution can be concatenated to get the desired solution. Refer to the following code:
% number of partitions of the time interval, tf .. final time
tspan = linspace(0,tf,N);
tAll = zeros(length(u_opt),1);
xAll = zeros(length(u_opt),4);
% Could be whatever the index of tspan you want
idx = 25;
for k = 1:numel(u_opt)
[tode,xode] = ode45(@(t,x)calc_f(t,x,u_opt(k),param),tspan,x0);
tAll(k) = tode(idx);
xAll(k,:) = xode(idx,:);