Can't stabilize system with PID

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:

[Sam Chak]

Hi @Andrei Rotaru

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]

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!

Answer 1

[Aquatris]

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)

Answer 2

[Sam Chak]

Hi @Andrei Rotaru

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

[Andrei Rotaru]

Thanks you! This is very helpful!