基于PID一阶倒立摆控制与仿真
目录
0. 前言 (1)
1 PID控制和编码器基本理论 (2)
2 方案设计 (3)
2.1 系统模型及分析 (3)
2.2 比例（P）控制 (4)
2.3 积分（I）控制 (4)
2.4 微分（D）控制 (4)
3 系统模型建立 (5)
3.1 一阶倒立摆的微分方程模型 (5)
3.2 传递函数 (7)
3.3 状态空间方程 (7)
4 软件编程与仿真 (8)
4.1 实际系统参数 (8)
4.2 PID控制设计分析 (8)
4.3 PID参数的确定 (10)
4.4 Simulink仿真 (11)
4.5 单极倒立摆建模 (13)
4.6 软件编程 (15)
5 系统调试和结果分析 (18)
6 结论及进一步设想 (19)
课设体会 (21)
基于PID一阶倒立摆控制与仿真
亲
摘要：本文主要研究目的是通过PID的调解实现对一阶倒立摆的控制，设计一个倒立摆的控制系统，使倒立摆这样一个不稳定的被控对象通过引入适当的控制策略使之成为一个能够满足各种性能指标的稳定系统。首先对平面一级倒立摆系统进行分析，根据具体参数建立数学模型，通过对模型的分析判断，设计倒立摆PID控制器。通过MATLAB软件进行仿真和实际系统实验，实现对倒立摆的稳定控制。但是由于PID控制器为单输入单输出系统，所以只能控制摆杆的角度，并不能控制小车的位置，所以小车会往一个方向运动。可以通过应用现代控制理论等单输入（小车加速度）多输出（摆杆角度和小车位置）的控制算法使系统更加的稳定。
关键词：倒立摆；PID控制；MATLAB仿真；
0. 前言
倒立摆是日常生活中许多重心在上、支点在下的控制问题的抽象模型，本身是一种自然不稳定体，它在控制过程中能有效地反映控制中许多抽象而关键的问题，如系统的非线性、可控性、鲁棒性等问题。对倒立摆系统的控制就是使小车以及摆杆尽快地达到预期的平衡位置，而且还要使它们不会有太强的振荡幅度、速度以及角速度，当倒立摆系统达到期望位置后，系统能克服一定范围的扰动而保持平衡。作为一种控制装置，它具有形象直观、结构简单、便于模拟实现多种不同控制方法的特点，作为一个被控对象它是一个高阶次、非线性、多变量、强耦合、不稳定的快速系统，只有采取行之有效的方法才能使它的稳定效果明了，因此对倒立摆的研究也成为控制理论中经久不衰的研究课题。倒立摆是进行控制理论研究的典型实验平台。倒立摆是机器人技术、控制理论、计算机控制等多个领域、多种技术的有机结合，其被控系统本身又是一个绝对不稳定、高阶次、多变量、强耦合的非线性系统，可以作为一个典型的控制对象对其进行研究。最初研究开始于二十世纪50 年代，麻省理工学院（MIT）的控制论专家根据火箭发射助推器原理设计出一级倒立摆实验设备。近年来，新的控制方法不断出现，人们试图通过倒立摆这样一个典型的控制对象，检验新的控制方法是否有较强的处理多变量、非线性和绝对不稳定系统的能力，从而从中找出最优秀的控制方法。
一阶倒立摆系统是一个典型的非线性、强耦合、多变量和不稳定系统。通过对它的研究不仅可以解决控制中的理论和技术实现问题，还能将控制理论涉及的主要基础学科：力学，数学和计算机科学进行有机的综合应用。其控制方法和思路无论对理论或实际的过程控制都有很好的启迪，是检验各种控制理论和方法的有效的“试金石”。倒立摆的研究不仅
有其深刻的理论意义，还有重要的工程背景。在多种控制理论与方法的研究与应用中，特别是在工程实践中，也存在一种可行性的实验问题，通过对倒立摆的控制，用来检验新的控制方法是否有较强的处理非线性和不稳定性问题的能力。同时，其控制方法在军工、航天、机器人和一般工业过程领域中都有着广泛的用途，如机器人行走过程中的平衡控制、火箭发射中的垂直度控制和卫星飞行中的姿态控制等
直线倒立摆是在直线运动模块上装有摆体组件，直线运动模块有一个自由度，小车可以沿导轨水平运动，在小车上装载不同的摆体组件，可以组成很多类别的倒立摆，直线柔性倒立摆和一般直线倒立摆的不同之处在于，柔性倒立摆有两个可以沿导轨滑动的小车，并且在主动小车和从动小车之间增加了一个弹簧，作为柔性关节。
1PID控制和编码器基本理论
在工程实际中，应用最为广泛的调节器控制规律为比例、积分、微分控制，简称PID控制，又称PID调节[6]。PID控制器问世至今已有近70年历史，它以其结构简单、稳定性好、工作可靠、调整方便而成为工业控制的主要技术之一。当被控对象的结构和参数不能完全掌握，或得不到精确的数学模型时，控制理论的其它技术难以采用时，系统控制器的结构和参数必须依靠经验和现场调试来确定，在积分控制中，控制器的输出与输入误差信号的积分成正比关系。对一个自动控制系统，如果在进入稳态后存在稳态误差，则称这个控制系统是有稳态误差的或简称有差系统。为了消除稳态误差，在控制器中必须引入“积分项”。积分项对误差取决于时间的积分，随着时间的增加，积分项会增大。这样即便误差很小，积分项也会随着时间的增加而加大，它推动控制器的输出增大使稳态误差进一步减小，直到等于零。因此，比例加积分(PI)控制器，可以使系统在进入稳态后无稳态误差。
倒立摆控制系统的工作原理是：由轴角编码器测得小车的位置和摆杆相对垂直方向的角度，作为系统的两个输出量被反馈至控制计算机。计算机根据一定的控制算法，计算出空置量，并转化为相应的电压信号提供给驱动电路，以驱动直流力矩电机的运动，从而通过牵引机构带动小车的移动来控制摆杆和保持平衡。绝对编码器通过与位数相对应的发光二极管和光敏二极管对输出的二进制码来检测旋转角度。与增量编码器原理相同，用于测量直线位移的传感器是光栅尺。由于光电编码器输出的检测信号是数字信号，因此可以直接进入计算机进行处理，不需放大和转换等过程，使用非常方便，因此应用越来越广泛。旋转编码器是一种角位移传感器，它分为光电式、接触式和电磁感应式三种，其中光电式脉冲编码器是闭环控制系统中最常用的位置传感器。
图1.1光电编码器原理示意图
旋转编码器有增量编码器和绝对编码器两种，图1.1为光电式增量编码器示意图，它由发光元件、光电码盘、光敏元件和信号处理电路组成。当码盘随工作轴一起转动时，光源透过光电码盘上的光栏板形成忽明忽暗的光信号，光敏元件把光信号转换成电信号，然后通过信号处理电路的整形、放大、分频、记数、译码后输出。为了测量出转向，使光栏板的两个狭缝比码盘两个狭缝距离小1/4节距，这样两个光敏元件的输出信号就相差π/2相位，将输出信号送入鉴向电路，即可判断码盘的旋转方向。
倒立摆系统的控制算法分两部分：摆杆自动起摆控制和稳定在垂直小角度范围内控制，一般以摆杆偏离垂直站立角度±12°为界限，在±12°范围内采用模糊控制算法使其稳定在垂直站立状态，超出这个范围则启动自动起摆控制算法，再次使其稳定在垂直站立状态。为了提高倒立摆系统的抗干扰性能，摆杆在±12°范围内摆动时使用的模糊规则库可进一步细分，以确定合适的控制电压。
2方案设计
2.1系统模型及分析
一阶倒立摆的模型如图2.1所示：
图2.1 一级倒立摆的模型示意图
倒立摆的物理构成可以表述为：光滑的导轨，可以在导轨上自由移动的小车，和一个质量块的摆杆。它们的铰接方式决定了它们在竖直平面内运动。水平方向的驱动力F使小车根据摆角的变化而在导轨上运动，从而达到倒立摆系统的平衡。根据设计要求，采用的方案如下。硬件部分由电机通过同步带驱动实现小车在滑杆上来回运动，保持摆杆平衡。电机编码器和摆角编码器向运动卡反馈小车和摆杆位置。
2.2比例（P）控制
比例控制是一种最简单的控制方式。其控制器的输出与输入误差信号成比例关系。当仅有比例控制时系统输出存在稳态误差（Steady-state error）。
2.3积分（I）控制
在积分控制中，控制器的输出与输入误差信号的积分成正比关系。对一个自动控制系统，如果在进入稳态后存在稳态误差，则称这个控制系统是有稳态误差的或简称有差系统（System with Steady-state Error）。为了消除稳态误差，在控制器中必须引入“积分项”。积分项对误差取决于时间的积分，随着时间的增加，积分项会增大。这样，即便误差很小，积分项也会随着时间的增加而加大，它推动控制器的输出增大使稳态误差进一步减小，直到等于零。因此，比例+积分（PI）控制器，可以使系统在进入稳态后无稳态误差。
2.4微分（D）控制
在微分控制中，控制器的输出与输入误差信号的微分（即误差的变化率）成正比关系。自动控制系统在克服误差的调节过程中可能会出现振荡甚至失稳。其原因是由于存在有较大惯性组件（环节）或有滞后（delay）组件，具有抑制误差的作用，其变化总是落后于误差的变化。解决的办法是使抑制误差的作用的变化“超前”，即在误差接近零时，抑制误差的作用就应该是零。这就是说，在控制器中仅引入“比例”项往往是不够的，比例项的作用仅是放大误差的幅值，而目前需要增加的是“微分项”，它能预测误差变化的趋势，这样，具有比例+微分的控制器，就能够提前使抑制误差的控制作用等于零，甚至为负值，从而避免了被控量的严重超调。所以对有较大惯性或滞后的被控对象，比例+微分（PD）控制器能改善系统在调节过程中的动态特性。
软件部分实现在 PID控制器各校正环节中，比例环节成比例地反映控制系统的偏差信号，偏差一旦产生，控制器立即产生控制作用，以减少偏差；积分环节主要用于消除稳态误差，提高系统的型别，积分作用的强弱取决于积分时间常数Ti，Ti越大，积分作用越弱，反之则越强；微分环节反映信号的变化趋势,即变化速率，并能在偏差信号值变得太大之前，在系统中引入一个有效地早起修正信号，从而加快系统的动作速度，减小调节时间。
PID控制器由比例单元(P)、积分单元(I)和微分单元(D)组成。其输入)(t e与输出)(t
u的
关系为：
ki
e
dt
kd
t
kpe
=(1.1)
+
u?+
t
de
t
*
)
(
/)
)
t
(
(dt
)
(
/1
)
式中积分的上下限分别是0和t，因此它的传递函数为：
比例作用下，通过现场试验找到等幅震荡的过渡过程，记下此时的比例度和等幅振荡周期，再通过简单的计算求出衰减振荡时控制器的参数。
3系统模型建立
3.1一阶倒立摆的微分方程模型
在忽略了空气流动，各种摩擦之后，可将倒立摆系统抽象成小车和匀质杆组成的系统，如下图2.2所示
图2.2 单级倒立摆模型示意图
下面我们对这个系统作一下受力分析。下图2.3是系统中小车和摆杆的受力分析图。其中，N和P为小车与摆杆相互作用力的水平和垂直方向的分量。
注意：在实际倒立摆系统中检测和执行装置的正负方向已经完全确定，因而矢量方向定义如图，图示方向为矢量正方向。
图2.3倒立摆模型受力分析
分析小车水平方向所受的合力，可以得到以下方程：
N x b F x M --=
（1） 由摆杆水平方向的受力进行分析可以得到下面等式：
)
sin (2
2θl x dt
d m
N +=
即 θθθθsin cos 2
ml ml x m N -+= （2） 把这个等式代入上式中，就得到系统的第一个运动方程：
F ml ml x b x
m M =-+++θθθθsin cos )(2 （3） 为了推出系统的第二个运动方程，我们对摆杆垂直方向上的合力进行分析，可以得到下面方程：
θθθθ
θcos sin )
cos (22
2 ml ml mg P l dt
d m
mg P --=-=-即： （4）
力矩平衡方程如下：
θ
θθ I Nl Pl =--cos sin （5） 注意：此方程中力矩的方向，由于θφθφφπθsin sin ,cos cos ,-=-=+=，故等式前面有负号。
合并这两个方程，约去P 和N ，由2
3
1ml I =得到第二个运动方程： θθθcos sin )(2x
ml mgl ml I -=++ （6） 设φπθ+=（φ是摆杆与垂直向上方向之间的夹角），假设φ与1（单位是弧度）相比很小，即φ《1，则可以进行近似处理：1cos -=θ，φθ-=sin ，0
)
(2
=dt
d θ。用u 来代表
被控对象的输入力F ，线性化后两个运动方程如下：
?????=-++=-+u
ml x b x m M x g ml I φφφ
)()(2 （7）
3.2 传递函数
对方程组(3)进行拉普拉斯变换，得到
?????=Φ-++=Φ-Φ+)
()()()()()()()()(2
22
22s U s s ml s s bX s s X m M s
s X s g s s ml I （8） 注意：推导传递函数时假设初始条件为0。
由于输出为角度φ，求解方程组（4）的第一个方程，可以得到
)(])
([)(2
2
s s
g l ml
ml I s X Φ-
+= （9）
角度与位置的传递函数：
mgl
s ml I mls
s X s -+=
Φ2
2
2)()
()( （10）
如果令x v
=，得到角度与速度的传递函数：
mgl
s ml I ml
s V s -+=
Φ2
2
)()
()( （11）
把上式代入方程组的第二个方程，得到
)()()()()()()(2
2222s U s s ml s s s g ml ml I b s s s g ml ml I m M =Φ-Φ?????
?+++Φ??????-++
（12） 整理后得到传递函数：
s
q
bmgl s
q
mgl
m M s q
bml s
s
q
ml
s U s -
+-
+
=Φ2
3
2
3
4
4
2
)()
()( （13）
其中 ])())([(22ml ml I m M q -++= （14） 3.3 状态空间方程
系统状态空间方程为
Du
CX y Bu AX X
+=+= （15）
方程组对φ
,x 解代数方程，得到解如下： ????
?
??
??++++++-==+++++-==u l m M l m M m M g x l m M b u
m M m M mg x m M b x x x )4(3)4()(3)4(3)4(4
)4(3)4(4φφ
φφφ （16）
整理后得到系统状态空间方程：
u l m M m M x x l
m M m M g l
m M b m M mg m M b x x
???
?????
????????+++?????
????????????????
??
??????
?+++-++-=
??
?
???????????)4(30)4(400)4()(3)4(3010
000)
4(3)
4(400010
φφφφ （17） u
x x x y ??????+??????
????????????
??=??????=0001
0000
01φφφ （18）
4 软件编程与仿真
4.1 实际系统参数
实际系统参数如下，求系统的传递函数、状态空间方程，并进行脉冲响应和阶跃响应的Matlab 仿真。
M 小车质量 1.096 Kg m 摆杆质量
0.109 Kg
b 小车摩擦系数
0 .1N/m/sec
l 摆杆转动轴心到杆质心的长度 0.2 5m
I 摆杆惯量
0.0034 kg*m*m
F 加在小车上的力
x 小车位置
θ 摆杆与垂直方向的夹角 T 采样频率
0.005秒
4.2 PID 控制设计分析
在一些情况下针对特定的系统设计的PID控制控制得很好，但它们仍存在一些问题需要解决：
如果自整定要以模型为基础，为了PID参数的重新整定在线寻找和保持好过程模型是较难的。闭环工作时，要求在过程中插入一个测试信号。这个方法会引起扰动，所以基于模型的PID参数自整定在工业应用不是太好。
如果自整定是基于控制律的，经常难以把由负载干扰引起的影响和过程动态特性变化引起的影响区分开来，因此受到干扰的影响控制器会产生超调，产生一个不必要的自适应转换。另外，由于基于控制律的系统没有成熟的稳定性分析方法，参数整定可靠与否存在很多问题。
因此，许多自身整定参数的PID控制经常工作在自动整定模式而不是连续的自身整定模式。自动整定通常是指根据开环状态确定的简单过程模型自动计算PID参数。
但仍不可否认PID也有其固有的缺点：
PID在控制非线性、时变、耦合及参数和结构不确定的复杂过程时，工作地不是太好。最重要的是，如果PID控制器不能控制复杂过程，无论怎么调参数都没用。
虽然有这些缺点，PID控制是最简单的有时却是最好的控制方法[8]
在计算机控制系统中，PID控制是通过计算机程序来实现的，因此它的灵活性很大。一些原来在模拟PID控制器中无法实现的问题，在引入计算机以后，就可以得到解决，于是产生了一系列的改进算法，形成非标准的控制算法，以改善系统的品质，满足不同的控制系统的需要。
这个控制问题和我们以前遇到的标准控制问题有些不同，在这里输出量为摆杆的位置，它的初始位置为垂直向上，我们给系统施加一个扰动，观察摆杆的响应。系统框图如下：
图2-2 直线一阶倒立摆闭环系统图
图中)
G是被控对象传递函数。
(s
(s
KD是控制器传递函数，)
考虑到输入0
r，结构图可以很容易地变换成
s
)
(
图2-3 直线一阶倒立摆闭环系统简化图
该系统的输出为
)
()
)(())(()
()()()
)(())((1)()
()(1)()(s F num numPID den denPID denPID num s y s F den denPID num numPID den num
s F s G s KD S G s y +=
+
=
+=
num ——被控对象传递函数的分子项 den ——被控对象传递函数的分母项 numPID ——控制器传递函数的分子项 denPID ——控制器传递函数的分母项
通过分析上式就可以得到系统的各项性能。由可以得到摆杆角度和小车加速度的传递函数：
mgl
s ml I ml
s V s -+=
Φ2
2
)()
()(
PID 控制器的传递函数为：
denPID
numPID s
K s K s K s K K K s KD I
P D I P D =
++=
+
+=2
)(
通过对控制量 v 双重积分即可以得到小车位置：
2
)()(s
s V s X =
4.3 PID 参数的确定
PID 控制器的参数整定是控制系统设计的核心内容。它是根据被控过程的特性确定PID 控制器的比例系数、积分时间和微分时间的大小。PID
控制器参数整定的方法很多，概括
起来有两大类：一是理论计算整定法。它主要是依据系统的数学模型，经过理论计算确定控制器参数。这种方法所得到的计算数据未必可以直接用，还必须通过工程实际进行调整和修改。二是工程整定方法，它主要依赖工程经验，直接在控制系统的试验中进行，且方法简单、易于掌握，在工程实际中被广泛采用。PID 控制器参数的工程整定方法，主要有临界比例法、反应曲线法和衰减法。三种方法各有其特点，其共同点都是通过试验，然后按照工程经验公式对控制器参数进行整定。但无论采用哪一种方法所得到的控制器参数，都需要在实际运行中进行最后调整与完善。现在一般采用的是临界比例法。利用该方法进行 PID 控制器参数的整定步骤如下：首先预选择一个足够短的采样周期让系统工作；其次仅加入比例控制环节，直到系统对输入的阶跃响应出现临界振荡，记下这时的比例放大系数和临界振荡周期；再次在一定的控制度下通过公式计算得到PID 控制器的参数。 4.4 Simulink 仿真
由实际系统的物理模型：
在 Simulink 中建立如图4所示的直线一级倒立摆模型：
图2.6 直线一级倒立摆PID 控制MATLAB 仿真模型
经过系统调试令Kp=200,Ki=10,K D=20，得到仿真结果如下：
26705
.00102125.002725
.0)
()
(2
-=Φs s V s
图2.7 PID控制参数
PID仿真结果：
图2.8 基于PID一阶倒立摆MATLAB仿真结果倒立摆摆杆角度曲线为：
图2.9 摆杆角度曲线
小车的位置输出曲线为：
图2.10 小车位置曲线
可以看出，由于PID 控制器为单输入单输出系统，所以只能控制摆杆的角度，并不能控制小车的位置，所以小车会往一个方向运动。 4.5 单极倒立摆建模
倒立摆系统的控制问题一直是控制研究中的一个经典问题。控制的目标是通过给小车底座施加一个力u （控制量），使小车停留在预定的位置，并使杆不倒下，即不超过一预先定义好的垂直偏离角度范围，图2.1为一级倒立摆的示意图，小车质量为M ，摆杆质量为m ，小车位置为x ，摆杆的角度为θ。
设摆杆偏离垂直线的角度为θ，同时规定摆杆重心坐标为),(.G G y x ，则 θs i n l x x G +=
θc o s
.l y G = 则根据牛顿定律，建立水平和垂直运动状态方程。 摆杆围绕其重心的转动运动可用力矩方程来描述
θθθ
cos sin Hl Vl I -=
摆杆重心的水平运动由下式描述
摆杆重心的垂直运动由下式描述
小车的水平运动有下式来描述
H
l x dt
d m =+)sin (2
2θmg
V l dt
d m -=θcos 2
2H
u x d M -=22
假设θ很小，sin θ≈θ，cos θ≈1.则以上各式变为
由上式可得出
由上式可得单极倒立摆方程
式中，
控制指标共有4个，即单极倒立摆的摆角θ、摆速θ 、小车位置x 和小车速度x 。将倒立摆运动方程转化为状态方程的形式。令x x x x x x ====)4(,)3(,)2(,)1(θθ，则状态方程为： Bu Ax x
+= 式中
对每个控制指标都选取PD 控制方式，控制器为
)()()(k dei k k ei k k u di PI I +=
式中，ei(k)和dei （k ）为控制指标i 的误差和误差变化率。
H u x
M mg V H l x m Hl Vl I -=-==+-= 0)(θθθ
θθθmgl x
ml ml l u ml x m M =++=++ )()(2
u
Mml
I m M ml
Mml
I m M gl M m m 2
2
)()()(++-
+++=θθ
u
Mml
l m M ml
I Mml
I m M gl
m x
2
2
2
2
2
)()(++++
++-=θ 2
1,12
12
=
=
l ml I ???????
??=00
10000000010
2
1
t t A ??????
?
??=43
00
t t B 2
2
42
32
2
222
1)(,)(,)(,)()(Mml
I m M ml
I t Mml
I m M ml
t Mml
I m M gl
m t Mml
I m M gl M m m t +++=
++-
=++=
+++=
∑==
4
1
)
()(l i
k u
k u
4.6软件编程
基于PID一阶倒立摆MATLAB仿真如下：clear all;
close all;
xiteP=100;
xiteI=30;
xiteD=30;
i1_1=-10;i2_1=-10;i3_1=-10;i4_1=-10;
p3_1=-10;p1_1=-10;p2_1=-10;p4_1=-10;
d1_1=-0;d2_1=-0;d3_1=-0;d4_1=-0;
error_1=0;
error_2=0;
e1_1=0;e2_1=0;e3_1=0;e4_1=0;
e1_2=0;e2_2=0;e3_2=0;e4_2=0;
u_1=0;
xk=[-10/57.3,0,0.020,0]; %初始状态
ts=0.02;
for k=1:1:1000
time(k)=k*ts;
Tspan=[0 ts];
para=u_1;
%para=0;
[t,x]=ode45('nn_pidf',Tspan,xk,[],para);
xk=x(length(x),:);
r1(k)=0.0; %摆角
r2(k)=0.0; %摆角速度
r3(k)=0.0; %汽车的位置
r4(k)=0.0; %汽车位率
x1(k)=xk(1);
x2(k)=xk(2);
x3(k)=xk(3);
x4(k)=xk(4);
e1(k)=r1(k)-x1(k);
e2(k)=r2(k)-x2(k);
e3(k)=r3(k)-x3(k);
e4(k)=r4(k)-x4(k);
error(k)=0.2*e1(k)+0.8*e3(k);
xx(1)=error(k)-error_1; %P
xx(2)=error(k); %I
xx(3)=error(k)-2*error_1+error_2; %D
i1(k)=i1_1+xiteI*xx(2);
p1(k)=p1_1+xiteP*xx(1);
d1(k)=d1_1+xiteD*xx(3);
i2(k)=i2_1+xiteI*xx(2);
p2(k)=p2_1+xiteP*xx(1);
d2(k)=d2_1+xiteD*xx(3);
i3(k)=i3_1+xiteI*xx(2);
p3(k)=p3_1+xiteP*xx(1);
d3(k)=d3_1+xiteD*xx(3);
i4(k)=i4_1+xiteI*xx(2);
p4(k)=p4_1+xiteP*xx(1);
d4(k)=d4_1+xiteD*xx(3);
de1(k)=e1(k)-e1_1;
ce1(k)=e1(k)-2*e1_1+e1_2;
u1(k)=i1(k)*e1(k)+p1(k)*de1(k)+d1(k)*ce1(k);
de2(k)=e2(k)-e2_1;
ce2(k)=e2(k)-2*e2_1+e2_2;
u2(k)=i2(k)*e2(k)+p2(k)*de2(k)+d2(k)*ce2(k);
de3(k)=e3(k)-e3_1;
ce3(k)=e3(k)-2*e3_1+e3_2;
u3(k)=i3(k)*e3(k)+p3(k)*de3(k)+d3(k)*ce3(k);
de4(k)=e4(k)-e4_1;
ce4(k)=e4(k)-2*e4_1+e4_2;
u4(k)=i4(k)*e4(k)+p4(k)*de4(k)+d4(k)*ce4(k);
%u(k)=u1(k)+u2(k)+u3(k)+u4(k);
u(k)=(11*u1(k)+20*u2(k)+96.5*u3(k)+35.5*u4(k))/50; if k==400
u(k)=u(k)+1
end
if u(k)>=10
u(k)=10;
elseif u(k)<=-10
u(k)=-10;
end
error_2=error_1;
error_1=error(k);
i1_1=i1(k);i2_1=i2(k);i3_1=i3(k);i4_1=i4(k);
p1_1=p1(k);p2_1=p2(k);p3_1=p3(k);p4_1=p4(k);
d1_1=d1(k);d2_1=d2(k);d3_1=d3(k);d4_1=d4(k);
e1_2=e1_1;e1_1=e1(k);
e2_2=e2_1;e2_1=e2(k);
e3_2=e3_1;e3_1=e3(k);
e4_2=e4_1;e4_1=e4(k);
u_1=u(k);
end
figure(1);
subplot(411);
plot(time,r1,'k',time,x1,'k'); %摆角
xlabel('time(s)');ylabel('Angle');
subplot(412);
plot(time,r2,'k',time,x2,'k'); %摆角速度
xlabel('time(s)');ylabel('Angle rate');
subplot(413);
plot(time,r3,'k',time,x3,'k'); %汽车位置 xlabel('time(s)');ylabel('Cart position'); subplot(414);
plot(time,r4,'k',time,x4,'k'); %汽车位率 xlabel('time(s)');ylabel('Cart rate'); figure(5);
plot(time,u,'k'); %力的变化 xlabel('time(s)');ylabel('Force'); nn_pidf.m 文件如下： function dx=dym(t,x,flag,para) global A B C D u=para; dx=zeros(4,1);
%State equation for one link inverted pendulum %dx=A*x+B*u;
dx=[x(2);(x(4)^2*sin(x(3))-3.67875*sin(2*x(3))-50.9784*x(2)+1)/(11-0.75*cos(x(3))^2)+para/(1+0.1*(1-0.75*cos(x(3))^2)); x(4);
(7.3575*sin(x(3))-0.75*cos(x(3))*(x(4)^2*sin(x(3))-3.67875*sin(2*x(3))-50.9784*x(2)+1)/(11-0.75*cos(x(3))^2))-0.75*cos(x(3))*para/(1+0.1*(1-0.75*cos(x(3))^2))];
5 系统调试和结果分析
运行后得到如下的仿真结果：
2468
101214161820
time(s)
F o r c e