有關代碼及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