倒立摆的LQR控制器算法的设计与仿真

1、LQR控制器算法的设计与仿真线性二次最优控制的概述应用经典控制理论设计控制系统，能够解决很多简单、确定系统的实际设计 问题。但是对于诸多新型而复杂的控制系统，例如多输入多输出系统与阶次较高 的系统，往往得不到满意的结果。这时就需要有在状态空间模型下建立的最优控 制策略。最优控制是现代控制理论的核心。 所谓最优控制，就是在一定条件下，在完 成所要求的控制任务时，使系统的某种性能指标具有最优值。根据系统不同的用 途，可提出各种不用的性能指标。最优控制的设计，就是选择最优控制，以使某 一种性能指标为最小。最优控制就是在满足一定的约束条件下， 寻求最优控制策略，使系统的某种 性能指标具有最优值的一种控

2、制方法。对于一个线性系统，如果其性能泛函是状 态变量或控制变量的二次型函数的积分，则这样的最优控制问题就称为线性二次 型最优控制问题，而利用线性二次型性能指标设计的控制器就被称作线性二次型 最优(Linear Quadratic Regulator)空制器，缩写为 LQR 控制器15。利用线性二次型最优控制算法不仅易于分析、处理和计算，得到的系统控制 方法还具有较好的鲁棒性与动态特性以及能够获得线性反馈结构等优点，而且特别可贵的是，线性二次型最优控制得到状态线性反馈的最优控制规律，易于构成闭环最优控制。而且，MATLAB的应用为线性二次型最优控制的理论仿真提供 了很好的条件，更为实现稳、准、快

3、的控制目标提供了方便。因而线性二次型控 制器已经成为自动控制系统中反馈控制设计的一种重要的控制方法，在实际的控制系统设计中得到了广泛的应用。线性二次最优控制的原理线性二次最优控制的控制原理框图如图：图4.1 LQR控制系统结构图Fig 4.1 Conventional LQR control system structure(4-1)对于状态方程式所表示的连续时间的线性被控对象有:x(t) Ax Bu(t),x(0)X0上式中，x(t)为n维状态向量；u(t)为m维控制向量；A, B分别为nxn及nxm 阶的常数矩阵。在进行线性二次最优控制系统设计时，比较令人感兴趣的是如何选择控制向量u(t)

4、,使得给定的性能指标达到最小，线性二次最优调节器(LQR)是针对系统状态方程，寻找最优控制，使得控制性能指标J达到最小，其中Q、 R分别表示了对状态变量和输入变量的加权值。二次型性能指标函数：(4-2)(4-3)(4-4)(4-5)(4-6)J (XTQX UTRU)dt则有如下状态反馈控制律：U(t) Kx(t)式中，K为最优反馈矩阵：K R1BTP在式(4-4)中,P为Riccati方程非负定对称解。而 Riccati方程为:_T _1_T_PA A P PBR B P Q 0如此，可得到状态反馈增益向量 K:K R1BTP由此可见最优控制器的设计的关键是选择合适的加权矩阵Q和R,并根据式

5、(4-6)可以算出P,这样就能求出反馈增益 K 了。而加权矩阵Q和R的具体作 用为：Q中某元素相对增加时，其对应的状态变量的响应速度增加，其他状态变量的响应速度相对减慢；R增加时，控制力减小，角度变化变小，跟随速度变慢。 改变矩阵Q的值，可以得到不同的响应效果，Q主对角线元素的值在一定范围 之内越大，系统调整时间越短，而且抵抗干扰的能力越强，但是Q不能过大，不然将对实验结果有一定的影响16o上述推导过程即为线性二次最优控制的控制原理。而当今现在，随着计算机 技术的飞速发展，已经可以不使用上述公式进行繁琐的计算，而利用 MATLAB 的lqr命令轻松的得到反馈矩阵K的值：K lqr(A,B,Q,

6、R)二次最优控制器的参数调整二次最优控制器的参数调整关键在于选择好合适的加权矩阵Q和R,这样就能得出反馈增益K17。其中，加权矩阵中的Q和R的具体作用为：Q中某元 素相对增加时，其对应的状态变量的响应速度增加， 其他状态变量的响应速度相 对减慢；R增加时，控制力减小，角度变化变小，跟随速度变慢。而在实际应用 中，通常将R值进行固定，然后对Q进行调整以得到较好的控制效果18。所以,在这里，选取R=0.003o在LQR控制器设计中，选取R=0.003, Q阵为对角阵:Q1,10000Q2,200Q00Q3,30000Q4,4其中，Q1,1代表小车位置的权重（即x）, Q2,2代表小车速度的权重(4

7、-7)(即 X)，Q3,3代表摆杆角度的权重（即），Q4,4代表摆杆角速度的权重（即）。而在该实验中，只考虑小车的位移及摆杆的摆动角度问题，所以，可以将 Q2,2和Q4,4都定义为0,这样，只需要调整Q1,1和Q3,3o则Q的矩阵可以变换为:Q1,10000000Q00Q3,300000(4-8)接下来，根据第二章所建立的系统的状态空间表达式:x 01x 00.002410000 x1.725470 x0101.00503 u000.00556 26.58199 02.32597(4-9)x1 0 0 0 x0 0 10(4-10)可以得到:0010.0024101.725470001.005

8、03AB0001000.0055626.5819902.3259710 0 0 C0 0 10那么，首先取Q1,1=Q3,3=1的时候，在MATLAB的输入界面输入以下程序: %最优控制A= 0 1 0 0;0 -0.00241 1.72547 00 0 0 1；0 -0.00556 26.58199 0;B= 0 1.00503 0 2.32597;C= 1 0 0 0;0 0 1 0；D= 0 0 ;%求向量KQ11=1;Q22=0;Q33=1;Q44=0;Q=Q11 0 0 0;0 Q22 0 00 0 Q33 0;0 0 0 Q44;R = 0.003;K = lqr(A,B,Q,R)

9、%计算LQR控制矩阵Ac = (A-B*K);Bc = B;Cc = C;Dc = D;%求阶跃响应XO=0.01 -0.001 0.001 -0.001; %初始状态t=0:0.005:5;U=20*ones(size(t)；Y,X=lsim(Ac,Bcn,Cc,Dc,U,t,XO);yl=Y(:,1);y2=Y(:,2);plot(t,yl,r-,t,y2,b-);legend(小车位置,摆杆角度)；axis(0 5 -30 30);grid可以得到如下的仿真结果：K =-18.26 -14.5166.3513.42系统的阶跃响应曲线如下图所示。图4.2 Q1,1=Q3,3=1时倒立摆系统

10、阶跃响应图Fig.4.2 When Q1,1 =Q3,3 = 1 inverted pendulum system step response figure由这个图可以看出，小车的位移和角度在 2.5秒的时候就基本能达到稳定， 而且就小车位移而言，开始时候产生的超调还是相对较小的。 但是，摆杆的角度, 在开始时候的超调量有些偏大。而倒立摆系统对摆杆角度的超调量的要求比较严 格，需要将摆杆角度的超调量控制在尽量小的范围内。又可根据Q值变化对系统的影响，增大Q3,3值可以减小超调量，但是同时也会增加系统的响应时间， 反之则超调量增大，响应时间减小；增大 Q1,1可以增大系统的超调量，而减小 系统的

11、响应时间，反之则超调量减小，响应时间增大。所以，在这里，为了达到 更好的实验效果，尝试增大 Q3,3值以减小系统的超调量，同时，考虑适当的增 大Q1,1来抵消Q3,3变化产生的响应时间增长的影响。经过多次尝试之后，在 Q1,1=5, Q3,3=80 时，得到：K=-47.46 -44.64 230.1234.07系统的阶跃响应曲线如下图所示：:片F-.r- i卜车位首-攫杆带度- n ajta t ! rLl/ :_ i-r - - ris _r 1X20to-ZD0 05 i 1 5 2 2 S 3 J 6 4 4 5 6图4.3 Q1,1=5, Q3,3=80时倒立摆系统阶跃响应图Fig.

12、4.3 When Q1,1= 5, Q3,3 = 80 inverted pendulum system step response figure由图(4.2)与图(4.3)对比，可以看出：当Q1,1和Q3,3增大时，系统摆杆角度的超 调量在减少，但是系统响应时间也会随之变长。因此，要使超调量较小，并具有 较短的响应时间，同时就要在可能的范围内，将Q1,1和Q3,3都选择的相对大一 点。因此，在本文中LQR控制器参数就选取 Q1,1=5, Q3,3= 80, R=0.003止匕时, 系统的的超调量和稳定时间都比较令人满意。直线一级倒立摆系统的LQR仿真根据直线一级倒立摆的状态空间方程模型，使用

13、MATLAB对倒立摆系统进行建模和仿真实验19。通过QUARC软件中的各种模块，建立了直线一级倒立 摆的最优控制仿真图如下图所示：图4.4直线一级倒立摆的LQR控制仿真模型Fig.4.4 LQR control simulation model of the linear inverted pendulum系统输入为小车的外力，状态变量是小车位移，摆杆角度，小车速度以及摆杆角速度这四个变量的实时信号，输出为小车位置 x,摆杆角度8信号，并将输出信号反馈给输入信号形成闭环。在本设计中，要求倒立摆系统保持竖直倒置的 摆杆的稳定，所以，系统的摆杆的初始角度，小车初始速度以及摆杆初始角速度 均为0。小车的位移输入是一个在第1秒产生阶跃的一个阶跃信号，系统通过阶 跃信号来使小车的位移发生变化，使小车左右移动，以此触发倒立摆系统的LQR 控制，控制摆杆的稳定。将系统的各项参数输入到仿真系统中，运行LQR控制器仿真系统。得出以下各结果：displacement 252015I 10 x 505012345678910time(s)图4.5小车位移仿真曲线angle1