Hello, kindly help. I am getting the following error: Error using odearguments @(T,X)STATEEQ(T,X,U,TU) must return a column vector. Error in Descent (line 17)

function

% Steepest descent method is used to find the solution.

eps = 1e-3;

options = odeset('RelTol', 1e-4, 'AbsTol',[1e-4]);

t0 = 0; tf = 1440;

% R = 0.1;

step = 0.4;

t_segment = 500;

Tu = linspace(t0, tf, t_segment); % discretize time

u = ones(1,t_segment); % guessed initial control u=1

initx = 0.1; % initial values for states

initp = 0; % initial values for costates

max_iteration = 100; % Maximum number of iterations

for i = 1:max_iteration

% 1) start with assumed control u and move forward

[Tx,X] = ode45(@(t,x) stateEq(t,x,u,Tu), [t0 tf], initx); %options);

% 2) Move backward to get the trajectory of costates

x = X(:,1);

[Tp,P] = ode45(@(t,p) costateEq(t,p,u,Tu,x,Tx), [tf t0], initp, options);

p = P(:,1);

% Important: costate is stored in reverse order. The dimension of

% costates may also different from dimension of states

% Use interploate to make sure x and p is aligned along the time axis

p = interp1(Tp,p,Tx);

% Calculate deltaH with x1(t), x2(t), p1(t), p2(t)

dH = pH(x,p,Tx,u,Tu);

H_Norm = dH'*dH;

% Calculate the cost function

E(i,1) = tf*((x')*x/length(Tx) + u/length(Tu));

% if dH/du < epslon, exit

if H_Norm < eps

% Display final cost

E(i,1)

break;

else

% adjust control for next iteration

u_old = u;

u = AdjControl(dH,Tx,u_old,Tu,step);

end;

end

% plot the state variables & cost for each iteration

figure(1);

plot(Tx, X ,'-');

hold on;

plot(Tu,u,'r:');

text(.2,0.08,'x(t)');

%text(.25,-.1,'x_2(t)');

text(.2,.4, 'u(t)');

s = strcat('Final cost is: E=',num2str(E(end,1)));

text(.4,1,s);

xlabel('time');

ylabel('states');

hold off;

print -djpeg90 -r300 eg2_descent.jpg

figure(2);

plot(E,'x-');

xlabel('Iteration number');

ylabel('J');

print -djpeg90 -r300 eg2_iteration.jpg

if i == max_iteration

disp('Stopped before required residual is obtained.');

end

% State equations

function dx = stateEq(t,x,u,Tu)

dx = zeros(2,1);

u = interp1(Tu,u,t); % Interploate the control at time t

R=50;

L=2;

q=10*sin(2*pi*Tu);

dx = 1/L*(q/u)-x*R/L-q/u;

% Costate equations

function dp = costateEq(t,p,u,Tu,x,xt)

dp = zeros(2,1);

x = interp1(xt,x,t); % Interploate the state varialbes

u = interp1(Tu,u,t); % Interploate the control

R=50;

L=2;

q=10*sin(2*pi*Tu);

dp = -2*x+p*R/L;

% Partial derivative of H with respect to u

function dH = pH(x,p,tx,u,Tu)

% interploate the control

u = interp1(Tu,u,tx);

R=50;

L=2;

q=10*sin(2*pi*Tu);

dH = 1+p*q/u^2-(p*q)/(L*u^2);

% Adjust the control

function u_new = AdjControl(pH,tx,u,tu,step)

% interploate dH/du

pH = interp1(tx,pH,tu);

u_new = u - step*pH;

Answer 1

[Cris LaPierre]

Your function stateEq returns a variable dx. The error is telling you it must be a colum vector. The most likely solution is to transpose dx.

dx = 1/L*(q/u)-x*R/L-q/u;

dx = dx';

[ILENIKEMANYA]

It doesn't solve the problem just yet.

[Cris LaPierre]

It fixes the error you are asking about. Please share more detais. Are you getting an error message? If so, please share the full error message (all the red text).

[ILENIKEMANYA]

yes, it fixed that one specifically.

I am now faced with these errors:

Error using odearguments

@(T,X)STATEEQ(T,X,U,TU) returns a vector of length 500, but the length of initial conditions vector is 1. The vector

returned by @(T,X)STATEEQ(T,X,U,TU) and the initial conditions vector must have the same number of elements.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in Descent (line 17)

[Tx,X] = ode45(@(t,x) stateEq(t,x,u,Tu), [t0 tf], initx, options);

[Cris LaPierre]

So you have a size mismatch. What is the intended size? 500 or 1? The output of your odefun should be the same size as your initial condition vector, which you named initx.

[Tx,X] = ode45(@(t,x) stateEq(t,x,u,Tu), [t0 tf], initx); %options);

[ILENIKEMANYA]

I need it to be 500, so how do I do that in that line?

[Cris LaPierre]

Then initx needs to a vector with 500 elements. It is currently a scalar. You need to define an initial condition for each state.

initx = 0.1;                    % initial values for states

If they all have the same initial condition, consider using ones to create the vector, and then multiply it by the initial condition.

[ILENIKEMANYA]

Why is this code just running nicely without a problem if their initial conditions are just two for each state and they have 50 elements?

% Edit to be able to call and run the function here
eg2OC_Descent

ans = 0.0282

function eg2OC_Descent
%EG2OC_Descent    Example 2 of optimal control tutorial.
%    This example is from D.E.Kirk's Optimal control theory: an introduction
%    example 6.2-2 on page 338-339.
%    Steepest descent method is used to find the solution.
eps = 1e-3;
options = odeset('RelTol', 1e-4, 'AbsTol',[1e-4 1e-4]);
t0 = 0; tf = 0.78;
% R = 0.1;
step = 0.4;
t_segment = 50;
Tu = linspace(t0, tf, t_segment);    % discretize time
u = ones(1,t_segment);               % guessed initial control  u=1
initx = [0.05 0];                    % initial values for states
initp = [0 0];                       % initial values for costates
max_iteration = 100;                 % Maximum number of iterations
for i = 1:max_iteration
    % 1) start with assumed control u and move forward
    [Tx,X] = ode45(@(t,x) stateEq(t,x,u,Tu), [t0 tf], initx, options);
    
    % 2) Move backward to get the trajectory of costates
    x1 = X(:,1); x2 = X(:,2);
    [Tp,P] = ode45(@(t,p) costateEq(t,p,u,Tu,x1,x2,Tx), [tf t0], initp, options);
    p1 = P(:,1);
    % Important: costate is stored in reverse order. The dimension of
    % costates may also different from dimension of states
    % Use interploate to make sure x and p is aligned along the time axis
    p1 = interp1(Tp,p1,Tx);
    
    % Calculate deltaH with x1(t), x2(t), p1(t), p2(t)
    dH = pH(x1,p1,Tx,u,Tu);
    H_Norm = dH'*dH;
    
    % Calculate the cost function
    J(i,1) = tf*(((x1')*x1 + (x2')*x2)/length(Tx) + 0.1*(u*u')/length(Tu));
    
    % if dH/du < epslon, exit
    if H_Norm < eps
        % Display final cost
        J(i,1)
        break;
    else
        % adjust control for next iteration
        u_old = u;
        u = AdjControl(dH,Tx,u_old,Tu,step);
    end; 
end
% plot the state variables & cost for each iteration
figure(1);
plot(Tx, X ,'-');
hold on;
plot(Tu,u,'r:');
text(.2,0.08,'x_1(t)');
text(.25,-.1,'x_2(t)');
text(.2,.4, 'u(t)');
s = strcat('Final cost is: J=',num2str(J(end,1)));
text(.4,1,s);
xlabel('time');
ylabel('states');
hold off;
print -djpeg90 -r300 eg2_descent.jpg
figure(2);
plot(J,'x-');
xlabel('Iteration number');
ylabel('J');
print -djpeg90 -r300 eg2_iteration.jpg
if i == max_iteration
    disp('Stopped before required residual is obtained.');
end
end

% State equations
function dx = stateEq(t,x,u,Tu)
dx = zeros(2,1);
u = interp1(Tu,u,t); % Interploate the control at time t
dx(1) = -2*(x(1) + 0.25) + (x(2) + 0.5)*exp(25*x(1)/(x(1)+2)) - (x(1) + 0.25).*u;
dx(2) = 0.5 - x(2) -(x(2) + 0.5)*exp(25*x(1)/(x(1)+2));
end
% Costate equations
function dp = costateEq(t,p,u,Tu,x1,x2,xt)
dp = zeros(2,1);
x1 = interp1(xt,x1,t);   % Interploate the state varialbes
x2 = interp1(xt,x2,t);
u = interp1(Tu,u,t);     % Interploate the control
dp(1) = p(1).*(u + exp((25*x1)/(x1 + 2)).*((25*x1)/(x1 + 2)^2 - ...
    25/(x1 + 2))*(x2 + 1/2) + 2) - ...
    2*x1 - p(2).*exp((25*x1)/(x1 + 2))*((25*x1)/(x1 + 2)^2 - ...
    25/(x1 + 2))*(x2 + 1/2);
dp(2) = p(2).*(exp((25*x1)/(x1 + 2)) + 1) - ...
    p(1).*exp((25*x1)/(x1 + 2)) - 2*x2;
end
% Partial derivative of H with respect to u
function dH = pH(x1,p1,tx,u,Tu)
% interploate the control
u = interp1(Tu,u,tx);
R = 0.1;
dH = 2*R*u - p1.*(x1 + 0.25);
end
% Adjust the control
function u_new = AdjControl(pH,tx,u,tu,step)
% interploate dH/du
pH = interp1(tx,pH,tu);
u_new = u - step*pH;
end

[ILENIKEMANYA]

Mine is just as it is there and It is running. That is the one I was trying to follow. My functions are all in one file.

[ILENIKEMANYA]

Now, what is the difference between that code and what I want to achieve? Now that It is running with 50 elements? though the initial condition is just one?

Kindly assist

[Cris LaPierre]

I've modified it so that it can run here. I suspect you can compare the differences between the two as easily as I can, so I'll point out a couple higher-level items.

In this code, the equations in the odefunc (e.g. stateEQ) return scalar values. In yours, it is returning arrays.

The code you have developed appears to be quite careless when it comes to the sizes of the variables used in calculations. This may be due to a misunderstaning of what ode45 is doing.

So this code is solving the equations at a specific interval in time, and ode45 walks through each time step. Your code, on the other hand, appears to be trying to solve the equations at all times at once, which is an incorrect approach when using an odesolver.

You might find value in going through some of the examples here: https://www.mathworks.com/help/matlab/ref/ode45.html

[ILENIKEMANYA]

I see what you are saying, thanks very much. It is caused by q=10*sin(2*pi*Tu).

q=10*sin(2*pi*t) is the correct one, my state is a function of time, I don't know how to treat that. My understanding is it will calculate q at each time. Or what am I suppose to do in that case? The other one was not a function of time.

Kindly advise again Please.

[Cris LaPierre]

You will know better than me if it is working correctly. Make the change, run the code, and inspect the results. Are they what you expect?