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

2 visualizzazioni (ultimi 30 giorni)
%% 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));
  1 Commento
Gonzalo Martínez de Pissón
Modificato: Gonzalo Martínez de Pissón il 9 Feb 2022
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!

Accedi per commentare.

Risposte (1)

Walter Roberson
Walter Roberson il 10 Feb 2022
%% 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
  2 Commenti
Rachel Ong
Rachel Ong il 10 Feb 2022
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
Walter Roberson il 10 Feb 2022
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) ?

Accedi per commentare.

Categorie

Scopri di più su Symbolic Math Toolbox in Help Center e File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by