bvp5c:Unable to solve the collocation equations -- a singular Jacobian encountered
8 views (last 30 days)
Show older comments
Hello MATLAB community,
I hope this message finds you well. I am relatively new to MATLAB and am currently working on solving a two-point boundary value problem using the bvp5c solver.
Unfortunately, I am facing challenges, and I would greatly appreciate any help or guidance.
Project Overview:
I am attempting to reprogram code for the optimal swing-up of a double inverted pendulum using the indirect method.
While the code in the paper "Optimal Swing up of Double Inverted Pendulum using Indirect Method" lacks MATLAB implementation details, I've made an attempt based on the project description.
My Code:
doublePendulumSwingUp()
function doublePendulumSwingUp
%parameters
J1 = 0.0024;
J2 = 0.0079;
l1 = 0.2;
l2 = 0.2;
a1 = 0.10;
a2 = 0.10;
m0 = 0.5;
m1 = 0.55;
m2 = 0.54;
c1 = 0.005;
c2 = 0.005;
g = 9.81;
% Set up time span
T = 2.2;
tspan = linspace(0, T, 1000);
initialConditions = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0];
solinit = bvpinit(tspan, initialConditions);
options = bvpset('RelTol', 1e-6, 'AbsTol', 1e-6);
sol = bvp5c(@odeFunction, @bcFunction, solinit, options);
figure;
plot(sol.x, sol.y);
title('Lösung des BVP-Problems');
xlabel('Time');
ylabel('State');
legend('w', 'phi1', 'phi2', 'dw', 'dphi1', 'dphi2', 'lam1', 'lam2', 'lam3', 'lam4', 'lam5', 'lam6');
function dydt = odeFunction(t, y)
w = y(1);
phi1 = y(2);
phi2 = y(3);
dw = y(4);
dphi1 = y(5);
dphi2 = y(6);
lam1= y(7);
lam2= y(8);
lam3= y(9);
lam4= y(10);
lam5= y(11);
lam6= y(12);
D11=m0+m1+m2;
D12=(m1*a1+m2*l1)*cos(phi1);
D13=m2*a2*cos(phi2);
D22=m1*a1^2+m2*l1^2+J1;
D23=m2*l1*a2*cos(phi1-phi2);
D33=m2*a2^2+J2;
b1=-(m1*a1+m2*l1)*dphi2^2*sin(phi1)-m2*a2*dphi2^2*sin(phi2);
b2=m2*l1*a2*dphi2^2*sin(phi1-phi2)+c1*dphi1-c2*(dphi2-dphi1);
b3=-m2*l1*a2*dphi1^2*sin(phi1-phi2)+c2*(dphi2-dphi1);
h1=0;
h2=(m1*a1+m2*l1)*g*sin(phi1);
h3=m2*a2*g*sin(phi2);
lam=[lam1;lam2;lam3;lam4;lam5;lam6];
u=-(lam4*(D22*D33-D23^2)-lam5*(D12*D33-D13*D23)-lam6*D13*D22-D12*D23)/(D11*D22*D33-D11*D23^2-D12^2*D33+2*D12*D13*D23-D13^2*D22);
u1=u;
u2=0;
u3=0;
M = [D11, D12, D13;
0, D22, D23;
0, 0, D33];
U=[u1;u2;u3];
b=[b1;b2;b3];
h=[h1;h2;h3];
xdot= M \ (U-b+h);
x1=w;
x2=phi1;
x3=phi2;
x4=dw;
x5=dphi1;
x6=dphi2;
x=[x1;x2;x3;x4;x5;x6];
%StateSpace Form
f1=x4;
f2=x5;
f3=x6;
f4=xdot(1);
f5=xdot(2);
f6=xdot(3);
f=[f1;f2;f3;f4;f5;f6];
H = 0.5*u + lam1*f1 + lam2*f2 + lam3*f3 + lam4*f4 + lam5*f5 + lam6*f6;
%Co State Equations
dlam1 = -gradient(H, w);
dlam2 = -gradient(H, phi1);
dlam3 = -gradient(H, phi2);
dlam4 = -gradient(H, dw);
dlam5 = -gradient(H, dphi1);
dlam6 = -gradient(H, dphi2);
dydt = [
f1;
f2;
f3;
f4;
f5;
f6;
dlam1;
dlam2;
dlam3;
dlam4;
dlam5;
dlam6];
end
function res = bcFunction(ya, yb)
res = [
ya(1);
ya(2) + pi;
ya(3) + pi;
ya(4);
ya(5);
ya(6);
yb(1);
yb(2);
yb(3);
yb(4);
yb(5);
yb(6)
];
end
end
3 Comments
Torsten
on 19 Jan 2024
Edited: Torsten
on 19 Jan 2024
As I said, in your computational flow, H is one single numerical value, not a function of the variables w, phi1, phi2, dw, dphi1 and dphi2.
Thus no operation (gradient,diff,...) can give you derivatives of H with respect to independent variables.
Maybe you could include your equations in a mathematical notation - it's easier to find a solution for your problem than reconstructing it from your code.
Answers (2)
Torsten
on 22 Jan 2024
I get
syms x1 x2 x3 x4 x5 x6 real
syms D11 D12 D13 D22 D23 D33 real
syms u l1 l2 l3 l4 l5 l6 f1 f2 f3 f4 f5 f6 h1 h2 h3 b1 b2 b3 f456 real
D = [D11 D12 D13;D12 D22 D23;D13 D23 D33];
f1 = x4;
f2 = x5;
f3 = x6;
f456 = inv(D)*([u;0;0]-[b1;b2;b3]-[h1;h2;h3]);
f4 = f456(1);
f5 = f456(2);
f6 = f456(3);
df4du = diff(f4,u);
df5du = diff(f5,u);
df6du = diff(f6,u);
solu = solve(u + l4*df4du + l5*df5du + l6*df6du==0,u);
f4 = subs(f4,u,solu);
f5 = subs(f5,u,solu);
f6 = subs(f6,u,solu);
u = solu;
H = 0.5*u^2+l1*f1+l2*f2+l3*f3+l4*f4+l5*f5+l6*f6;
dl1 = -diff(H,x1)
dl2 = -diff(H,x2)
dl3 = -diff(H,x3)
dl4 = -diff(H,x4)
dl5 = -diff(H,x5)
dl6 = -diff(H,x6)
dx1 = f1
dx2 = f2
dx3 = f3
dx4 = f4
dx5 = f5
dx6 = f6
18 Comments
See Also
Categories
Find more on Logical in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!