First of all, let's run a demo of inverted pendulum given in SIMULINK package.
- Preparation
Copy the model file for the demoC:\MATLAB\R2006b\toolbox\simulink\simdemos\simgeneral\penddemo.mdl
to your working directory, C:\MATLAB\R2006b\work.
- Execution of demo
On MATLAB command prompt, type as follows to open the model file of the demo:>>penddemo
Model window of penddemo
Start the simulation by select Simulation -> Start.
Animation window
In animation window, the triangle shows the reference position of the cart.
- Exercise:
- Move the horizontal bar to change the reference position.
- Open the signal generator
in the left hand side,
and change the values of Amplitude and/or Frequency. (default is zero.)
- Change mass of cart and/or length of pendulum as follows:
Open the block named `Pendulum'
by double-click, and specify
the parameters.
Is response changed? Does inverted pendulum fall down?
- Plot the response on figure
In the model file, there are blocks named `t' and `'y'.
Due to this blocks, time sequence is stored to variable t, and motion of the system (reference and actual position of cart, and pendulum angle) is stored.
In MATLAB prompt, type as follows to plot the graph:
>> plot(t, y)
Although the simulation period is set to infinity as default setting, you can stop when you want followed by execution of the plot command.
- Change gains of controller
The following part is the controller:
Controller
After double click the blocks, you can change its parameters.
- Kill controller
Just cut some signals and run simulation. Does pendulum fall down?
- Stop simulation
Select Simulation -> Stop to stop the simulation.
The equation of motion is defined in the block named `Pendulum'
.After single-clicking on Pendulum block, select Edit -> Look Under Mask to open the following window:
Contents of Pendulum block
In the above figure, the input U is force applied to the cart.
The cart position is outputted to the output port `Cart'.
The pendulum angle is outputted to the output port `Pendulum'
The equation of motions are described in two blocks named `f(u)'. Double click to open the contents.
The upper block of f(u) has input signals of force, angular velocity, and angle, and has output signals of acceleration velocity of cart..
The lower block of f(u) has output signal of acceleration velocity in angle.
In those blocks, force, angular velocity, and angle are shown as u(1), u(2), and u(3), respectively.
Exercise:
Now, check if the following equations of motion are described in two f(u) blocks:
For derivation procedure of the equation of motion, see this document (sorry now in Japanese) where the derivation procedure of equation of motion for 2-link inverted pendulum is shown. The case of 1-link can be derived by setting the mass of the second link to be zero.
The equation of motion in demo file above is derived by neglecting the moment of inertia of pendulum around its center of gravity.
By letting the moment of inertia to be J, the following equation of motion is derived.
Exercise:
Check if the equation of motion (3) is equivalent to (1) and (2) when J = 0.
The moment of inertia becomes as follow:
Exercise:
Modify the equation of motion described in two f(u) blocks so that the moment of inertia is considered.
Use variable J for moment of inertia.
Answer:
Fcn1 (acceleration of cart)
((m*l*l+J)*(u[1]+m*l*power(u[2],2)*sin(u[3]))-m*m*l*l*g*sin(u[3])*cos(u[3]))/((M+m)*(m*l*l+J)- power(m*l*cos(u[3]),2))Fcn2 (angular acceleration of pendulum)
((M+m)*m*g*l*sin(u[3])-(u[1]+m*l*power(u[2],2)*sin(u[3]))*m*l*cos(u[3]))/((M+m)*(m*l*l+J)-power(m*l*cos(u[3]),2))
Note:
変数 J を定義は，つぎのように行う:
pendulum ブロックを選択(左クリック)してから Edit タグ → Edit Mask → Initialization → Add と選択し，
Variable: フィールドに J を入力する．
Exercise
Compare the responses of two cases when J is given as equation (4) and when J is set to 0.
Example of execution:
(Double click the Reference block and specify the reference position
of cart as 5.)
(Then, run simulation with J = 0)
>>y0 = y;
(run simulation with J = 0.21*(0.61/2)^2/3)
>> plot(t, y0(:,2), 'b--', t, y(:,2), 'b', t, y0(:,3), 'r--', t, y(:,3), 'r')
>> xlabel('Time [sec]');
>> ylabel('Displacement of cart and angle of pendulum');
Red: pendulum angle
Broken line: J = 0
Line: J is given as equation (4)
Exercise:
Change the mass of cart and/or pendulum, and the length of pendulum, to compare the resultant responses before and after the changes.
Let's remove the all of controllers in demo file once, then re-construct our controller, which is state feedback controller.
In order to construct state feedback controller, first modify the Pendulum block so that it output the additional signals, which are velocity of cart and angular velocity of pendulum.
Next, remove controllers including PI controller, and create the structure of state feedback controller.
Use x - r as feedback signal of position, where x is the cart position and r is the reference position, in order for the cart would follow the reference position.
Control system with state feedback controller
Caution: the sign in summation point should be changed from the model file.
In addition, set F as the variable name in the Matrix Gain block.
So far, the structure of the control system has been created.
The remaining task is to determine F so that the closed loop system achieves desirable performance including stability.
Set the matrix F (in strictly speaking, F is a row vector with 4 scalars) by trial and error. It might be difficult very much to achieve stability.
Example (but not stable):
>> F = [1.5, 10, 1, 10];In order to change time step for animation, copy the m-file penddemo.mdl below to your directory. (This file is used to display animation by Animation block, since this file name is specified in the block.)
C:\Program Files\MATLAB\R2006b\toolbox\simulink\simdemos\simgeneral\pendan.mdlThen edit the line as below to change 0.1 to, for example, 0.01.
ts = [0.1 0];
In the below, you will design controller. Let the design specification is given as the following:
Design specification of controller:
- The closed-loop is stable. (of course!!)
- The settling time for cart position is less than 1 sec, where the settling time is the time after which the unit step response is closed within +/- 5% of the reference signal.
- The overshoot for cart position is less than 10%, where the overshoot is defined as (maximum value - reference value)/reference value * 100% in unit step response.
You can see equation(3) has non-linear term on theta.
Therefore, it is difficult to design controller and to analyze closed-loop stability.
This difficulties will be solved by linearization of the equation around equilibrium angle of theta. The result is as follows:
Moreover, (5) can be transformed into 1st differential equation instead 2nd one by extending the size of state vector:
where
You can easily check the closed-loop stability by using equation (6), sine the equation is linear on theta.
Actually, closed-loop is stable if and only if all eigenvalues of the matrix A + B F have negative real part, when state feedback controller f = F z is used.
Exercise:
Check if the closed-loop system via matrix F previously determined by trial and error, is stable or not.
Example of the answer: check.m
Example of execution:
>> check ans = -46.2033 0.3339 + 1.2311i 0.3339 - 1.2311i -0.9475In the above, there are two eigenvalues of which real part is positive, which means instability. (The response will be diverged with frequency 1.23 rad/sec corresponding the imaginary part of the unstable eigenvalues)
The matrix F can be calculated back from given desired closed-loop poles (eigenvalues). This technique is called pole placement.
For example, let us first consider that all of 4 poles are specified as -1. Then, corresponding matrix F is calculated as below:
Copy the following m-file to your working directory: m-file for pole placement: poleplace.m Example of execution:
>> F = poleplace(A, B, [-1, -1, -1, -1])
F =
0.0211 7.7639 0.0842 0.8598
>> eig(A + B*F)
ans =
-1.0002
-1.0000 + 0.0002i
-1.0000 - 0.0002i
-0.9998
As shown in the above, all 4 poles of closed-loop system are about 1,
and suitable gain matrix F is automatically determined.
Exercise:
Check if the closed-loop system has desired performance or not
m-file for performance evaluation: check_perf.m
Example of execution:
As you can see in the above figure, the closed-system is actually stable, however, the response is too slow.
What do we do in order to improve the response ?
- If we use -2 instead -1 ?
>> F = poleplace(A, B, [-2, -2, -2, -2])
- How about complex eigenvalues?
>> F = poleplace(A, B, [-1+i, -1-i, -1+i, -1-i])
Note:
In order for F to be determined as real matrix, all eigenvalues have to be complex conjugate.
Find suitable state-feedback gain F by comparing resultant time responses so that the given control specification would be satisfied.
Idea:
- The time response is mainly dependent to the eigenvalue whose absolute value of its real part is large.
- By introducing complex eigenvalues - a + b i, - a - b i (where a and b are positive real number), time response tends to have oscillation. As the result, overshoot might become larger, however, raising time might become smaller.
- When a is large, decay speed of response becomes fast.
- When b is large, frequency of oscillation becomes higher.
Example which satisfies the specification:
- All eigenvalues are set to same real number, and its
absolute value is set as large as possible.
>> F = poleplace(A, B, 6.1*[-1, -1, -1, -1])
Example which satisfies the specification
- Overshoot is utilized by introducing complex eigenvalues
>> F = poleplace(A, B, 3*[-1.5+i, -1.5-i, -1.5+i, -1.5-i])
Example which satisfies the specification
Note:
- When the absolute value of pendulum angle becomes large, the difference between actual plant and the linearized (approximated) model would become large. Therefore, the stabilizing controller for linearized plant might violate the closed-loop stability.
- It is impossible to apply too large force input to the cart since in the model file there is saturation block.