基於MATLAB的單級倒立擺模擬
有關程式碼及word文件請關注公眾號“浮光傾雲”，後臺回覆A010.02即可獲取
一、單級倒立擺概述
倒立擺是處於倒置不穩定狀態，人為控制使其處於動態平衡的一種擺，是一類典型的快速、多變數、非線性、強耦合、自然不穩定系統。由於在實際中存在很多類似的系統，因此對它的研究在理論上和方法上均有重要意義。
單級倒立擺系統（Simple Inverted Pendulum System）是由倒立擺和小車兩部分組成。小車依靠直流電動機施加控制力，可以在導軌上左右移動，其控制目標是在有限長導軌上使倒立擺能夠穩定豎立在小車上而不倒，達到動態平衡[1]
圖1 倒立擺小車系統
Fig1. Inverted Pendulum on a Cart
圖中F為施加於小車的水平方向的作用力， 是擺杆的傾斜角。若不給小車施加控制力，擺杆會左右傾斜，控制的過程即當擺杆出現偏角，在水平方向給小車以作用力，通過小車的水平運動，使擺杆保持在垂直位置，意即控制系統的狀態引數，以保持倒立擺的倒立穩定。
倒立擺系統由6大部分組成，並構成一個閉環結構，包括計算機、資料採集卡、電源及功率放大器、直流伺服電機、倒立擺本體和兩個光電編碼器等模組，如圖2所示。
圖2 倒立擺系統結構組成示意圖
Fig2 Structure of the Single Inverted Pendulum System
二、單級倒立擺數學模型
建立控制系統的數學模型有兩種基本方法，其一，機理建模法，即對系統各部分的運動機理進行分析，根據它們所依據的物理規律或化學規律分別列寫相應的運動方程，合在一起便成為描述整個系統的方程；其二，實驗建模方法，即人為地給系統施加某種測試訊號，記錄其輸出響應，並用適當的數學模型去逼近。主要用於系統運動機理複雜因而不便分析或不可能分析的情況[2]。
倒立擺的形狀較為規則，而且是一個絕對不穩定系統，無法通過測量頻率特性方法獲取其數學模型，故適合用數學工具進行理論推導。此處可以用牛頓力學法和拉格朗日方程法兩種方式進行推導。
2.1模型假設
為了建立倒立擺系統的數學模型，先作出如下假設：
（1） 倒立擺與擺杆均為勻質剛體；
（2） 各部分的摩擦力（力矩）與相對速度（角速度）成正比；
（3） 施加在小車上的驅動力與加在功率放大器上的輸入電壓成正比；
（4） 皮帶輪與傳送帶之間無滑動，傳送帶無伸長狀態；
（5） 訊號與力的傳遞無延時；
（6） 忽略空氣阻力。
2.2模型符號
（1）牛頓法模型符號
表1 牛頓力學法符號含義
|
符號 |
物理量 |
數值 |
|
M |
小車質量 |
0.5kg |
|
m |
擺杆質量 |
0.2kg |
|
l |
擺杆轉動軸心到擺杆質心的長度 |
0.3m |
|
I |
擺杆繞質心的轉動慣量 |
0.006 |
|
b |
小車滑動摩擦係數 |
0.1N/m/sec |
|
c |
擺杆轉動摩擦係數 |
0.1 |
|
F |
加在小車上的力 |
—— |
|
x |
小車位移 |
—— |
|
Φ |
擺杆與法向夾角 |
—— |
（2）拉格朗日方程法模型符號
表2 拉格朗日方程法符號含義
|
符號 |
物理量 |
數值 |
|
M |
小車質量 |
0.5kg |
|
m |
擺杆質量 |
0.2kg |
|
l |
擺杆轉動軸心到擺杆質心的長度 |
0.3m |
|
I |
擺杆繞質心的轉動慣量 |
0.006 |
|
b |
小車滑動摩擦係數 |
0.1N/m/sec |
|
c |
擺杆轉動摩擦係數 |
0.1 |
|
F |
加在小車上的力 |
—— |
|
x |
小車位移 |
—— |
|
擺杆與法線方向的夾角 |
—— |
|
|
q |
系統的廣義座標 |
—— |
|
T |
倒立擺系統的動能 |
—— |
|
D |
倒立擺系統的耗散能 |
—— |
|
V |
倒立擺系統的勢能 |
—— |
|
T0 |
倒立擺系統中小車的動能 |
—— |
|
D0 |
倒立擺系統中小車的耗散能 |
—— |
|
V0 |
倒立擺系統中小車的勢能 |
—— |
|
T1 |
倒立擺系統中擺杆的動能 |
—— |
|
D1 |
倒立擺系統中擺杆的耗散能 |
—— |
|
V1 |
倒立擺系統中擺杆的勢能 |
—— |
2.3數學模型
2.3.1 Newton-Euler法
對系統中的擺杆和小車分別進行受力分析，如圖3。其中，N和P為小車與擺杆相互作用力的水平和垂直方向的分量。 以順時針方向為向量方向，且角速度、角加速度都按同一方向為向量方向。
圖3 倒立擺模型受力分析
分析小車水平方向所受的力，可以得到以下方程：
(1)
由擺杆水平方向的受力進行分析，擺杆水平方向的合外力為N，則可以得到運動微分方程：
(2)
由式（2），得：
(3)
把式(3)代入式(1)中，得系統的第一個運動方程：
(4)
對擺杆垂直方向上的受力進行分析，可以得到下面方程：
(5)
即
(6)
以擺杆質心為基點，得擺杆的力矩平衡方程如下：
(7)
合併方程（3）、(6)和(7)，約去 P和N ，
得到系統的第二個運動方程：
（8）
進行近似處理，線性化後兩個運動方程如下：
（9）
因為擺體繞支點的轉動慣量I與擺體繞質心轉動慣量 關係為
（10）
將式（20）代入到式（19）中，得到：
（11）
2.3.2 Newton力學法
對系統中的擺杆和小車分別進行受力分析，如圖4。其中，N和P為小車與擺杆相互作用力的水平和垂直方向的分量。
圖4 倒立擺模型受力分析
分析小車水平方向所受的力，可以得到以下方程：
(1)
由擺杆水平方向的受力進行分析，擺杆水平方向需要保持平衡，故合外力為0，則可以得到平衡方程：
(2)
式中,
則
最後得第一個運動微分方程：
(4)
對擺杆垂直方向上的受力進行分析，可以得到下面方程：
(5)
式中，
得:
(6)
以擺杆質心為基點，得擺杆的力矩平衡方程如下：
(7)
合併方程（3）、(6)和(7)，約去 和 ，
得到系統的第二個運動方程：
（8）
則可以進行近似處理，線性化後兩個運動方程如下：
（9）
因為擺體繞支點的轉動慣量I與擺體繞質心轉動慣量 關係為
（10）
將式（20）代入到式（19）中，得到：
（11）
2.3.3 Lagrange方程法
Lagrange 方程是以能量觀點建立的運動方程式，需要從兩個方面分析，一個是表徵系統運動的動力學量——系統的動能和勢能，另一個是表徵主動力作用的動力學量——廣義力。
由n個關節部件組成的機械系統，其Lagrange方程應為：
（12）
其中，q為系統的廣義座標，表示系統中線位移和角度的變數；T為倒立擺系統的動能，V為倒立擺系統的勢能，D為倒立擺系統中的耗散能。
由圖4分析，可以得出，小車和擺杆的各部分能量表達式為：
圖5單級倒立擺分析圖
Fig 5. Simple Inverted Pendulum analysis chart
對於小車：
（13）
對於擺杆：
（14）
對於倒立擺系統，由式（13）和式（14），得：
， （15）
對小車而言，由式（15），得：
由式（12）則有：
（16）
而當 的時候，即對擺杆而言，由式（15），得：
由式（12）則有：
（17）
根據上面各式綜合，可以得到單級倒立擺系統動力學方程為：
（18）
進行線性化處理，根據分析可知，在接近平面位置時
因此，將式（18）簡化後得到動力學方程：
（19）
2.4模型分析
（1）微分方程
由上述分析，可以得到系統的微分方程為：
（11）
即
(20)
（2）傳遞函式
令初始條件為零，對於輸入為作用力F，輸出為擺杆角度時，對式（20）進行拉普拉斯變換，得到：
(21)
由於力F為輸入，角度 為輸出，求解以上方程組，得
（3）狀態空間方程
由現代控制理論，控制系統的狀態空間方程
令
則
故狀態空間方程為：
則有
三、單級倒立擺開環系統分析
3.1穩定性分析
（1）系統模型
由上節分析，系統的傳遞函式為：
利用MATLAB工具，求系統傳遞函式，程式程式碼見附錄1-CDHS.m，執行結果如下：
傳遞函式的多項式表示式為：
G =
0.06 s^2
-----------------------------------------------
0.0204 s^4 + 0.0724 s^3 - 0.4016 s^2 + 0.0588 s
由上節分析，可知系統的狀態空間方程為：
則有
利用MATLAB工具，求系統狀態方程，程式程式碼見附錄1-ZTFC.m，執行結果如下：
系統狀態方程各引數值為：
A =
x1 x2 x3 x4
x1 0 1 0 0
x2 0 -0.02837 0.417 -0.07092
x3 0 0 0 1
x4 0 -0.07092 25.54 -4.344
B =
u1
x1 0
x2 0.2837
x3 0
x4 0.7092
C =
x1 x2 x3 x4
y1 1 0 0 0
y2 0 0 0 0
y3 0 0 1 0
y4 0 0 0 0
D =
u1
y1 0
y2 0
y3 0
y4 0
（2）穩定性分析
由MATLAB程式碼分析系統的極點、並判斷穩定性，程式碼見附錄2-WDXFX.m：
系統極點：
0
-6.5986
2.8989
0.1507
有兩個極點位於右半平面：系統不穩定
3.2開環響應分析
（2）脈衝響應
MATLAB模擬的開環脈衝響應（即給系統加一個脈衝推力）曲線如圖6，程式程式碼見附錄2-MCXY.m
圖6 開環脈衝響應曲線
（3）階躍響應
MATLAB模擬的開環階躍響應曲線如圖7，程式程式碼見附錄2-JYXY.m，
圖7 開環階躍響應曲線
（5） Simulink模擬平臺分析
利用MATLAB中的simulink模擬平臺進行模擬，系統開環輸入為作用力F，輸出為擺杆角度，則將數值代入傳遞函式中，可得系統的傳遞函式值：
G =
0.06 s^2
-----------------------------------------------
0.0132 s^4 + 0.0724 s^3 - 0.4016 s^2 + 0.0588 s
在Simulink中搭建模擬平臺，如圖8所示：
圖8 Simulink模擬平臺
模擬執行，得到開環系統的階躍響應（a）和脈衝響應（b）如圖9：
（a）階躍響應曲線 （b）脈衝響應曲線
圖9 開環系統響應曲線
由開環響應得知，系統不穩定
四、單級倒立擺系統控制器設計及閉環系統分析
4.1 單級倒立擺系統的控制
倒立擺有多種控制方法，當前的控制方法可分為：
（1）線性理論控制方法
將倒立擺系統的非線性模型進行近似線性化處理，獲得系統在平衡點附近的線性化模型，然後再利用各種線性系統控制器設計方法，如PID控制、狀態反饋控制、LQR控制等，得到期望的控制器。
（2）非線性理論控制方法
由於線性控制理論與倒立擺系統多變數、非線性之間的矛盾，人們意識到針對多變數、非線性物件，採用具有非線性特性的多變數控制是解決多變數、非線性系統的必由之路。因此人們先後開展了倒立擺系統的非線性控制方法研究[3]。
4.2 PID控制器的設計
4.2.1 PID控制的基本原理
PID控制電路的主要原理是將偏差的比例（P）、積分（I）和微分（D）通過線性組合構成控制量，對被控物件進行控制。
圖10 PID控制框圖
PID控制器是一種線性控制器，它根據給定值與實際輸出值構成控制偏差
將偏差的比例（P）、積分（I）、微分（D）通過線性組合構成控制量，對被控物件進行控制，故稱為PID控制器，其控制規律為：
或寫成傳遞函式的形式：
簡單來說，PID控制器各校正環節作用如下[4]：
表3 各個控制環節的特點及作用
|
環節名稱 |
特點及作用 |
|
比例控制（P） |
減少了系統穩態誤差，提高了系統精度，但是降低了系統的相對穩定性。 |
|
積分控制（I） |
在系統前向通道加入積分環節，提高系統的型別，改善系統的穩定效能，但是，同樣會降低系統的相對穩定性，對系統的動態分析不利。 |
|
比例+微分控制（PD） |
比例加微分控制引入了開環零點，從而增加了系統的相角裕度，提高了系統的相對穩定性，改善了動態效能。 |
|
比例+積分+微分（PID） |
在比例+積分+微分的控制中，系統被引入了兩個零點，使系統的相角裕度大大提高，特別適用於大慣性環節。 |
4.2.2 引數設計
（1）系統分析
在開環系統中新增PID控制器構成閉環控制系統，系統框圖如圖11：
圖11 PID控制系統框圖
系統r=0，則可以將框圖變換為圖12，
圖12 變換後的系統框圖
構成閉環系統的傳遞函式為：
式中，賦值後得到傳遞函式的多項式表示式為：
G =
0.06 s^2
-----------------------------------------------
0.0204 s^4 + 0.0724 s^3 - 0.4016 s^2 + 0.0588 s
（2）引數選擇
篩選合適的PID引數，利用MATLAB中的PID Tuner模擬，當：引數為表4時達到較好的效能，相應的效能指標見表5
表4 較佳引數值
表5 效能指標
（3）程式碼分析
PID控制器傳遞函式為：
KD =
3.007 s^2 + 23.16 s + 36.42
---------------------------
s
構成的閉環系統傳遞函式為：
Gx =
0.1804 s^4 + 1.39 s^3 + 2.185 s^2
------------------------------------------------
0.0204 s^5 + 0.2528 s^4 + 0.9881 s^3 + 2.244 s^2
4.2.3閉環系統分析
（1）穩定性分析
MATLAB程式碼見附錄3-WDXFX.m，執行結果如下：
系統極點
0.0000 + 0.0000i
0.0000 + 0.0000i
-8.0836 + 0.0000i
-2.1543 + 2.9945i
-2.1543 - 2.9945i
系統穩定
（2）時域分析
用MATLAB程式碼（見附錄3—XNFX.m）分析系統的階躍響應，並求得效能指標如下：
終值C = 0.9737
峰值時間peak_time = 0.4558s
最大超調量max_overshoot = 51.1121%
上升時間rise_time = 0.1253s
調整時間settling_time = 1.9256s
階躍響應曲線如圖13
圖13 階躍響應曲線
（3）頻域分析
由MATLAB程式碼分析頻域特徵，繪製頻譜圖，如圖14所示，程式碼見附錄3-PYFX.m
圖14 頻域伯德圖
4.3 LQR控制器設計
現代控制理論的最突出特點就是將控制物件用狀態空間表示式的形式表示出來，這樣便於對多輸入多輸出系統進行分析和設計。線性二次型最優控制演算法（LQR）的目的是在一定的效能指標下，使系統的控制效果最佳，即利用最少的控制能量，來達到最小的狀態誤差，以下將利用最優控制演算法實現對一階倒立擺系統的擺杆角度和小車位置的同時控制。
此前的輸入是脈衝量，並且在設計控制器時，只對擺杆角度進行控制，而不考慮小車的位移。然而，對一個倒立擺系統來說，把它作為單輸出系統是不符合實際的，如果把系統當作多輸出系統的話，用狀態空間法分析要相對簡單一些。
4.3.1狀態空間系統分析
階躍響應曲線如圖15，分別是小車速度、位移和標杆角速度、角加速度，與外力F階躍輸入後的響應，程式程式碼見附錄4-ZTKJYXY.m
圖15 開環階躍響應曲線
可見，系統處於不穩定狀態。
4.3.2最優控制器的設計
最優控制的前提條件是系統是能控的，下面來判斷一下系統的能控能觀性。
圖16 狀態反饋原理圖
（1）線性系統的可控性
如果一個系統在有限的時間內，在某種控制向量u(t)的作用下，系統狀態向量x(t)可由一種初始狀態達到任意一種目標狀態時，則稱此係統是可控的。也就是說，一個系統是可控的，就是在理論上存在這樣的控制向量u(t)，使系統從當前狀態轉變為所希望的狀態。
規定初始時間
,初始狀態是任意的
,終端時間為
，終端狀態即目標狀態設在狀態空間的原點上，即 =0,可控性判據分為兩種形式：
方式一：設系統的狀態方程為：
式中，B暫時設為 列陣，u為輸入訊號
若 n×n矩陣
滿秩，即
則這樣的系統是可控的，否則不可控。
方式二：設系統的狀態方程為：
式中，B暫時設為 列陣，u為輸入訊號
如果方程式無重特徵根，那麼它可變成
式中， （i=1,2，···，n）是系統的相互不等的特徵根。
如果系統的輸入矩陣B中沒有一行是全為0時，則系統是可控的，否則不可控。
由MATLAB程式碼分析，知系統r=4，系統是可控的，程式見附錄4-KKKGFX.m。
Col =
0 0.2837 -0.0583 0.5173
0.2837 -0.0583 0.5173 -3.5483
0 0.7092 -3.1010 31.5899
0.7092 -3.1010 31.5899 -216.4684
r =
4
（2）線性系統的可觀測性
一般來說，用系統狀態變量表示的系統狀態不一定都能直接測出，而系統的輸出是必須能測出的變數。如何根據有限的輸出瞭解系統的狀態是系統可測性問題。系統的可測性的定義為：系統中如果每一個初始狀態 都可通過在一個有限時間間隔內由輸出量 的觀測量確定，則稱此係統是可觀測的。系統的可觀測性同樣有兩種方式判別：
方式一：設系統的狀態方程為：
y=Cx
式中，B暫時設為 列陣，u為輸入訊號
若矩陣
滿秩，即
則這樣的系統是可觀測的。
方式二：設系統的狀態方程為：
y=Cx
式中，B暫時設為 列陣，u為輸入訊號
如果方程式無重特徵根，那麼它可變成
y=Cx
式中， （i=1,2，···，n）是系統的相互不等的特徵根。
如果系統的輸入矩陣B中沒有一行是全為0時，則系統是可觀測的。
4.3.3 LQR控制器
用MATLAB中的lqr函式，可以得到最優控制器對應的k，要求用LQR控制演算法控制倒立擺擺動至豎直狀態，並可以控制倒立擺左移和右移。
G3 = -1.0000 -3.4148 89.8052 11.7563
小車受擾動後能恢復到平衡狀態，故符合要求，如圖
圖17 LQR控制器
五、模擬介面GUI動態模擬
GUI即圖形使用者介面（Graphical User Interface）的簡稱，又稱圖形使用者介面，是指採用圖形方式顯示的計算機操作使用者介面。GUI是借鑑並結合了電腦科學、美學、心理學、行為動作學及商業領域分析學的人機系統工程，將人—機—環境三者作為一個整體系統總體設計。與以往的計算機命令介面相比，GUI更容易讓使用者在視覺上得到接受，並獲得很好的使用者體驗。此次演示採用了GUI技術，實現了倒立擺控制的動畫演示。
5.1程式設計
採用最優LQR控制演算法來實現倒立擺小車的控制，主程式dlb.m包括以下幾個部分：
（1）模型引數的設定：採用 mc_CreateFcn()、mc_Ca33back()實現小車質量的設定，採用同樣的方式可以實現對擺杆長度和質量的設定。
（2）LQR引數設定：採用 _CreateFcn()和 _Callback()實現 ,i=4,3,2,1。採用 r_CreateFcn()和 r_Callback()實現 R。
（3）利用LQR計算K：由 lqrok_Callback(）完成。
（4）K 的設定：採用 _CreateFcn()和 _Callbak()實現 ,i=4,3,2,1。 。
（5）擺杆角度和小車位置初始值設定：實現小車水平位置設定與回撥 、小車水平拖動條設定與回撥，擺杆角度設定與回撥、擺杆角度拖動條的建立與回撥。
（6）干擾的輸入：干擾主要包括有衝擊、階躍和正弦等三種。
（7）模擬時間和步長的設定：主要利用Tedit__CreateFcn()和 Step__CreateFcn()來實現。
（8）模擬啟動的設定。
（9）重置按鈕：利用 reset_Callback(),實現倒立擺模型引數和 LQR 引數的重置。
（10）退出按鈕：exit_callback()。
5.2演示介面的設計
首先在 MATLAB 環境命令框內輸入“guide”便可進入 GUI 設計介面。此次設計的 GUI 介面檔案為 dlb.fig，建立並儲存該檔案的同時，系統會自動生成主程式框架 dlb.m。 在主程式 dlb.m 框架下，通過相應的 GUI 組建回撥函式中描述模型和編寫控制演算法。通過 MATLAB 環境下執行“guideDLB.fig”即可開啟 GUI 編輯介面。 此次演示設計中，採用了 GUI 開發環境的觸控按鈕、靜態文字、可編輯文字、滑動條、 座標軸等系統元件。以小車的水平位置的介面設計為例，首先建立小車水平位置的 Edittext， 將其屬性 tag 標籤定義“in_po”。基於 GUI 的倒立擺 LQR 控制演示主介面如圖 18 所示。
圖18 GUI演示介面
5.3演示過程
通過以下4個步驟即可實現對倒立擺的動畫演示：
（1）輸入倒立擺的引數：包括倒立擺擺杆質量 m 和長度 L，還有小車質量 M。
（2）通過“初始值”下的“水平位置”和“擺杆角度”可設定倒立擺和小車的初始角度及位置；選擇一種干擾輸入（正弦、衝擊、階躍），並設定干擾輸入的幅值。
（3）根據控制系統要求的效能輸入控制器設計引數 和 R，點選“確定”即可得到控制器增益“K”。
（4）點選“啟動模擬”便可以實現倒立擺的動態模擬演示。
5.4模擬結果
（1）根據之前的計算結果取 M=0.5，m=0.2，l=0.3 即下圖 L=0.6，隨意設定小車初始時刻的水平位置和擺杆的角度。控制器的引數設定為 和R=1，利用 LQR 方法求長，計算可得控制器增益
（2）K=[-92.3048 -11.5154 -10 -13.2052].
（3）設定模擬時間為 10 秒，步長 0.1。
（4）點選“啟動模擬”圖示，可得到如下模擬結果如圖 19 所示。
圖19 模擬結果顯示
六、結論
在此次模擬中，對單級倒立擺小車系統分別採用牛頓力學法和拉格朗日方程法進行數學建模，並得到了相應的狀態方程，二者相比，拉格朗日方程從能量的角度進行建模更加簡便，便於處理。
對單級倒立擺小車系統的控制分別採用經典控制理論設計控制器（PID控制器）和現代控制理論設計控制器（最優控制LQR控制器），並用 MATLAB R2017b 軟體模擬和除錯控制系統。用MATLAB對不同控制方式下的倒立擺系統進行模擬，最後利用圖形介面GUI工具設計倒立擺的模擬效果動畫演示。小車位置跟蹤輸 入訊號，控制擺杆超調足夠小，穩態誤差滿足要求，上升時間和穩定時間也符合應用指標。
參考文獻
[1]郝麗娜,劉興剛.計算機模擬技術及CAD.[M].高等教育出版社.2009.
[2]呂慶莉.單級倒立擺系統分析與設計.陝西科技大學學報.2009,27(6):121-128.
[3]李洪興,苗志巨集,王加銀.四級倒立擺系統的變論域自適應模糊控制.中國科學E輯.2002,32(1):65-75.
[4]黃國盛.單級倒立擺小車的控制系統.[M].2017.09