Multiple PID tuning in order to control all four states in the inverted pendulum model

Hi!

I have to control the four states of the classical nonlinear inverted pendulum on a cart model (position and velocity of the cart, angle and angular velocity of the pendulum) in Simulink through PID control. Being on a cart it doesn't have to go through the swing up, the initial condition for the angle is π-0.5 (desired angle π with a small perturbation), so it's just the balancing problem.

This is my Simulink control scheme:

My implementation works as it should, but I had to manually tune the four PIDs because I wasn't able to obtain the same system with just one controller (if I give a vector of the four errors as input to a single PID block designed with a vector of four gains it messes up and creates a 4-dimensional output instead of executing the row-column product) and I tried everything but couldn't tune the multiple PIDs at the same time. Of course, the tuner app embedded in every PID block is useless in my case since the output of the system isn't just depending on a single PID's control input. I wonder if I can automatically find the optimal choices for the PID controllers, can someone help me?

[Sam Chak]

Hi @Claudio Muzi,

To really get a hand on this, we'll need to see what's going on inside the invPend() function. Without having access to the underlying model, designing a meaningful PID controller is going to be tricky. In that case, you'd pretty much have to resort to some good old-fashioned trial-and-error tuning to get it dialed in.

[Claudio]

Hi @Sam Chak, this is the code of the InvPend model:

If needed, the initial condition is , the desired states are and the gains I have manually found are the following:

[Sam Chak]

Hi @Claudio Muzi

Thanks for showing the "Code". But I cannot copy the code from the image and test it.

[Claudio]

Sorry for the misunderstanding, here's the code of the model:

function eq_dot = InvPend(x,v,theta,omega,u)
% Parameters
m=1; % mass of the pendulum
M=5; % mass of the cart
L=2; % length of the pendulum
g=-9.80665; % gravitational acceleration
d=1; % friction coefficient

% Equations
D=m*L^2*(M+m*(1-cos(theta)^2));
x_dot=v;
v_dot=(1/D)*(-m^2*L^2*g*cos(theta)*sin(theta)+m*L^2*(m*L*omega^2*sin(theta)-d*v))+m*L^2*u/D;
theta_dot=omega;
omega_dot=(1/D)*((m+M)*m*g*L*sin(theta)-m*L*cos(theta)*(m*L*omega^2*sin(theta)-d*v))-m*L*cos(theta)*u/D;
eq_dot=[x_dot; v_dot; theta_dot; omega_dot];
end

Answer 1

[Sam Chak]

Hi @Claudio Muzi

Given that the initial angle is less than from the equilibrium inverted position, the sine of the angle () can be considered a relatively weak nonlinearity within that range. As a result, the inverted pendulum can be treated as a linear system and a full-state feedback control scheme can be confidently applied in this case:

.

The control gains can be tuned through a trial-and-error approach. For this simple demonstration, integer values were used.

x   = pi-0.5:0.0001:pi+0.5;
y   = sin(x);
plot(x, [pi - x; y]), grid on
legend('Linear', 'Nonlinear')
title('Nonlinearity of sin(\theta)')
xlabel('Angle / radian'), ylabel('sin(\theta)')

Simulation in MATLAB

%% Inverted Pendulum on a Cart
function EoM = InvPend(t, x)
    % parameters
    m   =  1;       % mass of the pendulum
    M   =  5;       % mass of the cart
    L   =  2;       % length of the pendulum
    g   = -9.80665; % gravitational acceleration
    d   =  1;       % friction coefficient
    D   =  m*L^2*(M + m*(1 - cos(x(3))^2));

    % full-state feedback controller
    K   = [-1 -6 152 64];   % control gains
    u   = - K*[x(1);
               x(2);
               x(3) - pi;
               x(4)];

    % Equations of Motion
    Dx      = x(2);
    Dv      = (1/D)*(- m^2*L^2*g*cos(x(3))*sin(x(3)) + m*L^2*(m*L*x(4)^2*sin(x(3)) - d*x(2))) + m*L^2*u/D;
    Dtheta  = x(4);
    Domega  = (1/D)*((m + M)*m*g*L*sin(x(3)) - m*L*cos(x(3))*(m*L*x(4)^2*sin(x(3)) - d*x(2))) - m*L*cos(x(3))*u/D;
    EoM     = [Dx; 
               Dv; 
               Dtheta;
               Domega];
end

tspan   = [0 30];
x0      = [0; 0; pi-0.5; 0];
[t, x]  = ode45(@InvPend, tspan, x0);
plot(t, x), grid on, xlabel('t'), title('Inverted Pendulum on a Cart')
legend({'$x$', '$v$', '$\theta$', '$\omega$'}, 'Interpreter', 'LaTeX', 'fontsize', 14)

%% Steady-state value of theta
ss_theta = x(end,3)

ss_theta = 3.1416

The pendulum's angular position (orange curve) is maintained at the desired equilibrium at steady-state.

[Claudio]

First of all, thank you for the detailed answer, much appreciated.

Yeah, the linearization works, but what if I must deal with the nonlinear model? Can the control system tuner app be useful in this case? If the answer is yes, how? Otherwise, are there other ways?

[Sam Chak]

Hi @Claudio Muzi

I did not mention the linearization work in my previous response. If you look at the code (specifically the Equations of Motion section), you'll see that I dealt with your original nonlinear dynamics directly, but is stabilized using a full-state feedback controller.

In fact, this approach is equivalent to utilizing four Proportional-only Controllers. You can maintain the existing PID blocks in your Simulink model, but adjust the values in the "P" fields accordingly.

The Control System Tuner may be a useful tool. This app can linearize the pendulum model at the desired operating point, allowing you to tune linear controllers using graphical and automated methods. However, it maybe not user-friendly if the user is unfamiliar with most of the control terminology in the app. For example, the app expects user to know and specify the characteristics of the desired response model. I hope the app allow non-control users to directly specify "what-I-wish-is-what-I-get" parameters in future release.

Anyhow, I would recommend following the tutorial for further guidance.

https://www.mathworks.com/videos/automatic-tuning-of-a-helicopter-flight-control-system-90590.html

[Claudio]

Ok, so basically I can't use any tool similar to the embedded PID block tuner because I have to stick to the four PIDs (or just the proportional part, but still four components) representation?

Plus, I watched the tutorial but it seems to me I'm missing something. The probability it is just me not using the app in the right way and misinterpreting something is very high, but I can't even manage to get the manually found gains. Could you show me how to do it?

[Sam Chak]

Hi @Claudio Muzi

The PID block tuner is built for handling Single-Input-Single-Output (SISO) systems. However, the inverted pendulum on a cart is an underactuated Single-Input-Multiple-Output (SIMO) system, which means only a single controller is required.

In my previous comment, I was simply stating that the single controller is equivalent to putting four Proportional values in the corresponding PID blocks .

However, the Autotuner feature in the PID block is likely not helpful in this design.

When dealing with nonlinear MIMO systems, I usually advise writing out the mathematical formulation of the controller and then designing it in MATLAB using whatever custom tools you have learned from programming books or practical nonlinear control experts. If successful, only then should you implement the controller in Simulink.

[Claudio]

Perfect, got it. What about the control system tuner instead? Could it be used to tune the four PIDs? I watched the tutorial you posted but I can't really adapt it to my system

[Sam Chak]

I seldom use the Control System Tuner. That's why I posted the tutorial. As mentioned above, the gains are tuned through a trial-and-error, which I prefer to Search & Sort the gains (within a certain range) until the desired response is found. I believe Computer Scientists call this "brute-force search".

[Claudio]

The trial-and-error works fine, but it's not the optimal solution nor necessarily something pretty close to it. My question is: is it possible to determine the best (or a very close one) gains combination in terms of overshoot, settling time and control input cost? I thought the control system tuner could do something like that

Answer 2

[Sam Chak]

Hi @Claudio Muzi

Your new question in this comment is related to optimal control. I'm not sure if this topic is covered in your class (it wasn't in mine). However, it is possible to determine the optimal gains subject to a custom design cost function that you wish to optimize. One approach is to mentally linearize the nonlinear system (which I did below - quite straightforward with explanations), and then find the gains using the Linear Quadratic Regulator (LQR) algorithm.

Why LQR? Because it has a pre-designed, fixed cost function that is convex (quadratic) in nature. This means the control designers don't need to expend significant effort applying complicated mathematics to design a convex cost function.

Code to apply lqr() command to find the optimal gains (took me nearly two hours to search LQR info, links and type this out)

%% Parameters
m   =  1;       % mass of the pendulum
M   =  5;       % mass of the cart
L   =  2;       % length of the pendulum
g   = -9.80665; % gravitational acceleration
d   =  1;       % friction coefficient

%% Linear State-Space
Dpi = m*L^2*(M + m*(1 - cos(pi)^2));    % Comes from "D" and theta is set to pi

% These 4 elements are determined from x_dot (1st state equation)
a11 = 0;                            % zero because no x
a12 = 1;                            % coefficient of v
a13 = 0;                            % zero because no theta
a14 = 0;                            % zero because no omega

% These 4 elements are determined from v_dot (2nd state equation)
a21 = 0;                            % zero because no x
a22 = - (m*L^2*d)/Dpi;              % coefficient of v
a23 = - (m^2*L^2*g*cos(pi)^2)/Dpi;  % coefficient of sin(theta), d/dθ sin(θ) = cos(θ)
a24 = 0;                            % zero because the term is multiplied by omega = 0

% These 4 elements are determined from theta_dot (3rd state equation)
a31 = 0;                            % zero because no x
a32 = 0;                            % zero because no v
a33 = 0;                            % zero because no theta
a34 = 1;                            % coefficient of omega

% These 4 elements are determined from omega_dot (4th state equation)
a41 = 0;                            % zero because no x
a42 = (m*L*cos(pi)*d)/Dpi;          % coefficient of v
a43 = (m + M)*m*g*L*cos(pi)/Dpi;    % coefficient of sin(theta), d/dθ sin(θ) = cos(θ)
a44 = 0;                            % zero because the term is multiplied by omega = 0

% State matrix
A   = [a11 a12 a13 a14 
       a21 a22 a23 a24
       a31 a32 a33 a34
       a41 a42 a43 a44]

A = 4x4

0 1.0000 0 0 0 -0.2000 1.9613 0 0 0 0 1.0000 0 -0.1000 5.8840 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

b1  = 0;                            % zero because no u in x_dot
b2  = (m*L^2)/Dpi;                  % coefficient  of u in v_dot
b3  = 0;                            % zero because no u in theta_dot
b4  = - m*L*cos(pi)/Dpi;            % coefficient  of u in omega_dot

% Input matrix
B   = [b1
       b2
       b3
       b4]

B = 4x1

0 0.2000 0 0.1000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% Linear Quadratic Regulator
Q   = eye(size(A));
R   = 1;
K   = lqr(A, B, Q, R)

K = 1x4

-1.0000 -5.5777 151.7848 63.9705

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

The optimal gains determined by LQR are nearly the same as the ones I brute-force searched by trial-and-error using integers. However, I'm uncertain how to obtain these LQR values using the Simulink-based Control System Tuner app.

If the MATLAB-based solution I've provided is still unacceptable, I would suggest consulting experts on the Control System Tuner tool.

[Claudio]

Thanks a lot for the time and effort you put in your answers. I am familiar with LQR, but I don't know how to compute it in Simulink, and unfortunately in this case I have to use Simulink and PID controllers. I read about algorithm such as Bees Algorithm or Genetic Algorithm that have been successfully used in the inverted pendulum PID control, but I really don't know anything about it.

By the way, how and where can I contact someone who can use well the Control System Tuner tool?

[Sam Chak]

Hi @Claudio

OP: I am familiar with LQR, but I don't know how to compute it in Simulink and unfortunately in this case I have to use Simulink and PID controllers.

Simulink is a technically sophisticated graphical interface for modeling and simulating dynamical systems; however, it cannot solve optimal control problems on its own. In other words, you will need interactive tools from Control System Toolbox, Simulink Control Design, and Simulink Design Optimization to collaborate effectively for optimization-based tuning of controller gains.

Original Question: I wonder if I can automatically find the optimal choices for the (four) PID controllers, can someone help me?

Since you are familiar with LQR, you should be able to compute the gains without relying on specific computational software. My proposed solution (2nd Answer) is to obtain a linear state-space model and then use LQR to find the optimal gains. This can be accomplished in MATLAB, Mathematica, Maple, Python, etc. Once the gains are obtained, these values can be implemented in the desired control architecture as, .

OP: I read about algorithm such as Bees Algorithm or Genetic Algorithm that have been successfully used in the inverted pendulum PID control, but I really don't know anything about it.

Implementing Artificial Bee Colony (ABC) or Genetic Algorithms (GA) involves "searching" for gains in MATLAB and sending the computed gains back to Simulink for running the simulation. The simulation results are then fed back to MATLAB for analysis. If the optimal value of the cost function is found, the "search" stops. Otherwise, the ABC or GA algorithms in MATLAB continue searching until a termination condition is met.

If your goal is strictly to achieve this process solely within Simulink and its associated tuning tools (such as the PID Autotuner), then the use of ABC and GA methods is not feasible.

OP: Is it possible to determine the best (or a very close one) gains combination in terms of overshoot, settling time and control input cost?

I think I understand what you're looking for. You're interested in using the Automated Control Design feature in Simulink to meet specific performance requirements, similar to what's demonstrated in this example:

https://www.mathworks.com/help/sldo/gs/optimize-controller-parameters-to-meet-step-response-requirements-gui.html

However, the challenge you're encountering is tuning all four PID Controller blocks simultaneously.

OP: where can I contact someone who can use well the Control System Tuner tool?

You might want to ask @Paul. I learned a lot from him about miscellaneous linear control design topics when I first started around the time of Covid.

[Paul]

I'm not at all familiar with the Control System Tuner tool.

[Claudio]

At this point, I'm satisfied with one way of autotuning the PIDs which gives me four triplets of gains

[Sam Chak]

Hi @Claudio,

Before carrying out the tuning task, you should check if the following 4-PID control architecture is desirable:

and

.

[Claudio]

If you mean acceptable, yes it is.

[Sam Chak]

Since the 4-PID control architecture is desirable, your next step would to simplifying the 12-parameter control equation such that , , , and so on. If similar terms are grouped, you can cut the number of parameters by half, and it will generally improves the tuning efficiency of your chosen optimization strategy.

[Oren]

The original question asked about PID control, but this answer only provides the proportional gain. How does one find the integral and differential gain? The PID Tuner app provided in matlab reccomends splitting the 4 states into 4 separate signals for 4 different PID blocks and tuning them separately. But tuning just assesses its response to a unit step function, and the velocity states don't converge to one (because that would be unstable), so the tuner app concludes they are broken.