PID model for MIMO system

Hi! Does anyone know how to make a PID controller (Matlab, Simulink) for the model from this article: https://ethz.ch/content/dam/ethz/special-interest/mavt/dynamic-systems-n-control/idsc-dam/Research_DAndrea/Cubli/Cubli_IROS2012.pdf

Answer 1

[Sam Chak]

Hi @Andrei Rotaru

Since the Cubli system is underactuated and you seek to use only a PID Controller, you must select which output you would like to control. This is necessary because PID control design generally works for single-input, single-output (SISO) systems.

Once a stabilizing PID controller is successfully designed, it can then be implemented on the multi-input, multi-output (MIMO) system. However, there is no guarantee that the other two outputs will be stable as well. Nevertheless, since the second output depends on the first output, if the first output is stable, the second output will also be stable.

%% 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   = [1 0 0];
sys = ss(A, B, C, 0*C*B);

%% Transfer function model
Gp  = tf(sys)

Gp =
 
        -5.214 s + 5.117e-17
  --------------------------------
  s^3 + 0.31 s^2 - 99.36 s - 8.717
 
Continuous-time transfer function.

%% PID Controller
Gc  = pidtune(Gp, 'PIDF')

Gc =
 
             1            s    
  Kp + Ki * --- + Kd * --------
             s          Tf*s+1 

  with Kp = -38.1, Ki = -51.7, Kd = -3.95, Tf = 0.00051
 
Continuous-time PIDF controller in parallel form.

%% Closed-loop transfer function
Gcl = minreal(feedback(Gc*Gp, 1))

Gcl =
 
           4.06e04 s^2 + 3.896e05 s + 5.284e05
  -----------------------------------------------------
  s^4 + 1960 s^3 + 4.111e04 s^2 + 1.949e05 s + 5.113e05
 
Continuous-time transfer function.

sys = ss(Gcl);
op  = findop(sys, y=1)

op = 
  OperatingPoint with properties:

            x: [4x1 double]
            u: 0.9677
            w: [0x1 double]
           dx: [4x1 double]
            y: 1
           rx: [4x1 double]
           ry: 0
    Equations: 5
     Unknowns: 5
       Status: 'Well-posed problem. Successfully computed the unique solution.'

%% Response to Custom Step Input
opt             = RespConfig;
opt.InputOffset = 0;
opt.Amplitude   = op.u;
step(Gcl, opt), grid on

[Andrei Rotaru]

Thank you! This is great!