Anyone knows how to solve this size compatibility issue? I got this from an example and I am trying to understand the script.

%% Solution Using Symbolic Toolbox (STB) in MATLAB Version 6.0
%% Example 2.14 pg 77 of NAIDU optimal control book --> not working 

clear all
S = dsolve('Dx1=x2, Dx2=-lam2, Dlam1=0, Dlam2=-lam1, x1(0)=1, x2(0)=2, x1(tf)=3, lam2(tf)=0')
t = 'tf'

eq1 = subs(S.x1)-'x1tf'
eq2 = subs(S.x2)-'x2tf'
eq3 = subs(S.lam1)-'lam1tf'
eq4 = subs(S.lam2)-'lam2tf'
eq5 = str2sym('lam1tf*x2tf-0.5*lam2tf^2')
S2 = solve(eq1, eq2, eq3, eq4, eq5)

% 'tf,x1tf,x2tf,lam1tf,lam2tf','lam1tf<>0'
%% lam1tf<>0 means lam1tf is not equal to 0;
%% This is a condition derived from eq5.
%% Otherwise, without this condition in the above
%% SOLVE routine, we get two values for tf (1 and 3 in this case)
tf = S2.tf
x1tf = S2.x2tf
x2tf = S2.x2tf
clear t
x1 = subs(S.x1)
x2 = subs(S.x2)
lam1 = subs(S.lam1)
lam2 = subs(S.lam2)

%% Convert the symbolic values to
%% numerical values as shown below.
j = 1;
tf = double(subs(S2.tf))

%% ERROR IS FROM HERE ONWARDS 

%% Converts tf from symbolic to numerical
for tp = 0:0.05:tf
    t = sym(tp);
    x1p(j) = double(subs(S.x1))
    x2p(j) = double(subs(S.x2))
    up(j) = -double(subs(S.lam2))
    tl(j)=tp;
    j = J+1;
end

plot(tl, x1p, tl, x2p, tl, up)
xlabel('t')
legend('x_1(t)', 'x_2(t)', 'u(t)')

The error says:

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in ex2_14 (line 41)

x1p(j)=double(subs(S.x1));

[Gonzalo Martínez de Pissón]

Hello!

It seems that the output of subs(S.x1) it's a vector of 2 elements and you are saving in one element. That's why the size mismatch happens. Try initializing each vector (x1p, x2p, up) and concatenating its solution as it is shown below:

%initialize x1p, x2p and up
x1p=[];x2p=[]; up=[];
for tp=0:0.05:tf
t=sym(tp);
%% coverts tp from numerical to symbolic
x1p=[x1p double(subs(S.x1))];
%% subs substitutes S.xl to xlp
x2p=[x2p double(subs(S.x2))];
%% double converts symbolic to numeric
up=[up -double(subs(S.lam2))];
%% optimal control u = -lambda_2
tl(j)=tp;
j=j+1 ;
end

hope it works!

Answer 1

[Walter Roberson]

%% Solution Using Symbolic Toolbox (STB) in MATLAB Version 6.0
%% Example 2.14 pg 77 of NAIDU optimal control book --> not working 

syms x1(t) x2(t) lam1(t) lam2(t)
syms l tf
Dx1 = diff(x1); Dx2 = diff(x2); Dlam1 = diff(lam1); Dlam2 = diff(lam2);
eqn = [Dx1==x2, Dx2==-lam2, Dlam1==0, Dlam2==-lam1, x1(0)==l ,x2(0)==2] % [x1(tf)==3, lam2(tf)==0]

eqn(t) = 

S = dsolve(eqn)

S = struct with fields:
      x2: (C1*t^2)/2 + C2*t + 2
      x1: (C1*t^3)/6 + (C2*t^2)/2 + 2*t + l
    lam1: C1
    lam2: - C2 - C1*t

x1vars = setdiff(symvar(S.x1), [t, l])

x1vars = 

bc = subs([S.x1, S.lam2], t, tf)

bc = 

bcsol = solve(bc == [3, 0], x1vars(1:2))

bcsol = struct with fields:
    C1: (3*(l + 2*tf - 3))/tf^3
    C2: -(3*(l + 2*tf - 3))/tf^2

Sbc = subs(S, bcsol)

Sbc = struct with fields:
      x2: (3*t^2*(l + 2*tf - 3))/(2*tf^3) - (3*t*(l + 2*tf - 3))/tf^2 + 2
      x1: l + 2*t - (3*t^2*(l + 2*tf - 3))/(2*tf^2) + (t^3*(l + 2*tf - 3))/(2*tf^3)
    lam1: (3*(l + 2*tf - 3))/tf^3
    lam2: (3*(l + 2*tf - 3))/tf^2 - (3*t*(l + 2*tf - 3))/tf^3

syms x1tf x2tf lam1tf lam2tf
assume(lam1tf ~= 0)
eq1 = Sbc.x1 - x1tf

eq1 = 

eq2 = Sbc.x2 - x2tf

eq2 = 

eq3 = Sbc.lam1 - lam1tf

eq3 = 

eq4 = Sbc.lam2 - lam2tf

eq4 = 

eq5 = lam1tf*x2tf-0.5*lam2tf^2;
S2 = solve([eq1,eq2,eq3,eq4,eq5])

S2 = struct with fields:
       l: [2×1 sym]
       t: [2×1 sym]
      tf: [2×1 sym]
    x1tf: [2×1 sym]
    x2tf: [2×1 sym]

% tf,x1tf,x2tf,lam1tf,lam2tf','lam1tf<>0')

%% lam1tf<>0 means lam1tf is not equal to 0;
%% This is a condition derived from eq5.
%% Otherwise, without this condition in the above
%% SOLVE routine, we get two values for tf (1 and 3 in this case)
%%
tf = S2.tf

tf = 

x1tf = S2.x1tf;
x2tf = S2.x2tf;
x1 = subs(Sbc.x1, S2)

x1 = 

x2 = subs(Sbc.x2, S2)

x2 = 

lam1 = subs(Sbc.lam1, S2)

lam1 = 

lam2 = subs(Sbc.lam2, S2)

lam2 = 

%% Convert the symbolic values to
%% numerical values as shown below.

%% coverts tf from symbolic to numerical
tf = S2.tf

tf = 

But that is as far as you can get. tf cannot be converted into an integer without knowing lam1tf, which you have not defined.

tpvals = 0:0.05:double(max(tf));

Error using : (line 38)
Unable to compute number of steps from 0 to -2*(1/lam1tf)^(1/2) by 1/20.

num_tp = length(tpvals)
for j = 1 : num_tp
    tp = tpvals(j);
    x1p(j, :) = double(subs(x1, t, tp));
    x2p(j, :) = double(subs(x2, t, tp));
    up(j, :) = -double(subs(lam2, t, tp));
    tl(j) = tp;
end

[Rachel Ong]

Why is it that now that the structure is changed everything does not have numbers? Previoulsy, I could still get some numbers for lam1tf just that the back part could not be subbed in properly

[Walter Roberson]

S = dsolve('Dx1=x2, Dx2=-lam2, Dlam1=0, Dlam2=-lam1, x1(0)=1, x2(0)=2, x1(tf)=3, lam2(tf)=0')

Warning: Support of character vectors and strings will be removed in a future release. Use sym objects to define differential equations instead.

S = struct with fields:
      x2: (3*t^2*(tf - 1))/tf^3 - (6*t*(tf - 1))/tf^2 + 2
      x1: 2*t - (3*t^2*(tf - 1))/tf^2 + (t^3*(tf - 1))/tf^3 + 1
    lam1: (6*(tf - 1))/tf^3
    lam2: (6*(tf - 1))/tf^2 - (6*t*(tf - 1))/tf^3

That does not define lam1tf

eq3 = subs(S.lam1)-'lam1tf'
eq5 = str2sym('lam1tf*x2tf-0.5*lam2tf^2')

Those use lam1tf, but do not define it.

% 'tf,x1tf,x2tf,lam1tf,lam2tf','lam1tf<>0'
%% lam1tf<>0 means lam1tf is not equal to 0;
%% This is a condition derived from eq5.
%% Otherwise, without this condition in the above
%% SOLVE routine, we get two values for tf (1 and 3 in this case)

Comments that mention lam1tf, nothing executable there

... and nothing else even mentions it, so nothing can solve for it or define it.

Perhaps you meant to define specific variables to solve() the list of equations for, instead of relying on MATLAB to choose 5 variables (for 5 equations) ?