Can't stabilize system with PID

8 views (last 30 days)
Hi. I'm trying to stabilize a model for this articol: https://ethz.ch/content/dam/ethz/special-interest/mavt/dynamic-systems-n-control/idsc-dam/Research_DAndrea/Cubli/Cubli_IROS2012.pdf , with a PID controller, and it's not working. I don't even know that this is the corect order of the cascade PID. I even try to take each individual part of the cascade PID and simulate it and use the tuner for that. And still it doesn't working. Does anyone know what to do?
This is the cascade PID:
This is the attempt to try to take individual part of the cascade and simulate:
This is the Control System tuner:
  2 Comments
Sam Chak
Sam Chak on 26 Jun 2024
Edited: Sam Chak on 26 Jun 2024
First, you should show us how you obtained the transfer function in the image. Using parameters given in the paper (code provided by @Aquatris), you can compare and notice that there are discrepancies. Most likely, you have supplied your own values to the parameters.
%% Parameters
l = 0.085;
lb = 0.075;
mb = 0.419;
mw = 0.204;
Ib = 3.34e-3;
Iw = 0.57e-3;
Cb = 1.02e-3;
Cw = 0.05e-3;
Km = 25.1e-3;
g = 9.81;
%% State-space model
A = [0 1 0
(mb*lb + mw*l)*g/(Ib + mw*l^2), -Cb/(Ib + mw*l^2), Cw/(Ib + mw*l^2)
-(mb*lb + mw*l)*g/(Ib + mw*l^2), Cb/(Ib + mw*l^2), -Cw*(Ib + Iw + mw*l^2)/(Iw*(Ib + mw*l^2))];
B = [0
- Km/(Ib + mw*l^2)
Km*(Ib + Iw + mw*l^2)/(Iw*(Ib + mw*l^2))];
C = eye(3);
sys = ss(A, B, C, 0);
%% Transfer functions
G = tf(sys)
G = From input to output... -5.214 s - 1.791e-16 1: -------------------------------- s^3 + 0.31 s^2 - 99.36 s - 8.717 -5.214 s^2 - 1.922e-16 s + 1.658e-18 2: ------------------------------------ s^3 + 0.31 s^2 - 99.36 s - 8.717 49.25 s^2 + 9.33 s - 4376 3: -------------------------------- s^3 + 0.31 s^2 - 99.36 s - 8.717 Continuous-time transfer function.
%% For SISO systems only (may not work on Cubli)
Gp = G(1);
Gc = pidtune(Gp, 'pidf')
Gc = 1 s Kp + Ki * --- + Kd * -------- s Tf*s+1 with Kp = -36.9, Ki = -49, Kd = -3.91, Tf = 0.000521 Continuous-time PIDF controller in parallel form.
Gcl = feedback(Gc*Gp, 1);
step(Gcl), grid on
Andrei Rotaru
Andrei Rotaru on 27 Jun 2024
Hi! Yees, I change the values with my own. It was pretty hard to understand what to do, but after your comment and some other research it's working. Thanks a lot!

Sign in to comment.

Accepted Answer

Aquatris
Aquatris on 26 Jun 2024
Edited: Aquatris on 26 Jun 2024
One thing you can do is use symbolic PID gains, plug them in your closed loop transfer function, and find the values via pole placement strategy. Below I provide the code where I use full state feedback and find a proportional gain to stabilize it, kind of like an LQR without optimization.
% model Parameters
l = 0.085;
lb = 0.075;
mb = 0.419;
mw = 0.204;
Ib = 3.34e-3;
Iw = 0.57e-3;
Cb = 1.02e-3;
Cw = 0.05e-3;
Km = 25.1e-3;
g = 9.81;
% model state-space matrices
A = [0 1 0;
(mb*lb+mw*l)*g/(Ib+mw*l^2) -Cb/(Ib+mw*l^2) Cw/(Ib+mw*l^2);
-(mb*lb+mw*l)*g/(Ib+mw*l^2) Cb/(Ib+mw*l^2) -Cw*(Ib+Iw+mw*l^2)/(Iw*(Ib+mw*l^2))];
B = [0
-Km/(Ib+mw*l^2)
Km*(Ib+Iw+mw*l^2)/(Iw*(Ib+mw*l^2))];
C = eye(3); % all outputs [theta_b theta_bDot theta_wDot]
sys = ss(A,B,C,0); %state space representation
% symbolic gain K, which is 1x3 matrix, full state feedback
syms K [1 3]
closedLoopPoles = eig(A-B*K*C); % for stability, we want [ real(closedLoopPoles) < 0 ]
% Gains found from pole placement, poles should be at -8 -20 and -27
sol = solve(closedLoopPoles==[-8;-20;-27],K);
Warning: Possibly spurious solutions.
% 2nd elements of the solution are actual numbers so use them
Kval = double([sol.K1(2),sol.K2(2),sol.K3(2)])
Kval = 1x3
-196.5043 -19.8323 -0.9892
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
% form the closed loop and check stability
closedLoop = feedback(sys*Kval,eye(3));
damp(closedLoop)
Pole Damping Frequency Time Constant (rad/seconds) (seconds) -8.00e+00 1.00e+00 8.00e+00 1.25e-01 -2.00e+01 1.00e+00 2.00e+01 5.00e-02 -2.70e+01 1.00e+00 2.70e+01 3.70e-02
t = 0:0.001:10;
u = [ones(length(t),1) zeros(length(t),1) zeros(length(t),1)]; % inputs for [theta_b theta_bDot theta_wDot]
lsim(closedLoop,u,t)

More Answers (1)

Sam Chak
Sam Chak on 27 Jun 2024
Here is the alternative control solution for finding the gains using LQR approach.
%% Parameters
l = 0.085;
lb = 0.075;
mb = 0.419;
mw = 0.204;
Ib = 3.34e-3;
Iw = 0.57e-3;
Cb = 1.02e-3;
Cw = 0.05e-3;
Km = 25.1e-3;
g = 9.81;
%% State-space model
A = [0 1 0
(mb*lb + mw*l)*g/(Ib + mw*l^2), -Cb/(Ib + mw*l^2), Cw/(Ib + mw*l^2)
-(mb*lb + mw*l)*g/(Ib + mw*l^2), Cb/(Ib + mw*l^2), -Cw*(Ib + Iw + mw*l^2)/(Iw*(Ib + mw*l^2))];
B = [0
- Km/(Ib + mw*l^2)
Km*(Ib + Iw + mw*l^2)/(Iw*(Ib + mw*l^2))];
C = eye(3);
sys = ss(A, B, C, 0);
%% LQR gain design
K = lqr(A, B, eye(size(A)), 1)
K = 1x3
-216.4821 -22.5414 -1.0020
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
%% Closed-loop system
closedLoop = feedback(sys*K, eye(size(A)))
closedLoop = A = x1 x2 x3 x1 0 1 0 x2 -1029 -117.7 -5.214 x3 1.056e+04 1110 49.25 B = u1 u2 u3 x1 0 0 0 x2 1129 117.5 5.224 x3 -1.066e+04 -1110 -49.35 C = x1 x2 x3 y1 1 0 0 y2 0 1 0 y3 0 0 1 D = u1 u2 u3 y1 0 0 0 y2 0 0 0 y3 0 0 0 Continuous-time state-space model.
%% Simulation
t = 0:0.01:2;
u = zeros(length(t), 3);
x0 = [pi/4, 0, 0]';
lsim(closedLoop, u, t, x0), grid on

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!