LQR控制系统优化

引言


在Habr上发表了几篇文章[1,2,3],这些文章直接或间接地与此主题相关。 在这方面,大家一定会注意到标题为“手指上的数学:线性二次调节器”的出版物[1],该出版物通俗地解释了最佳LQR控制器的工作原理。

在研究了动态优化方法的实际应用之后,我想继续这个主题,但是已经在使用Python的一个具体示例中进行了探讨。 首先,简要介绍一下动态优化的术语和方法。

优化方法分为静态和动态。 控制对象在各种外部和内部因素的影响下处于连续运动的状态。 因此,对控制时间T给出控制结果的评估,这是一个动态优化任务。

使用动态优化的方法,解决了在一定时间内分配有限资源的问题,并将目标函数写为一​​个积分函数。

解决此类问题的数学仪器是变分方法:经典的变分演算,最大L.S.原理。 Pontryagin和动态编程R. Bellman。

在时间,频率和状态空间中执行控制系统的分析和综合。 状态空间中控制系统的分析和综合被引入课程中,但是,使用SQR控制器的培训材料中介绍的技术是为Matlab设计的,不包含分析的实际示例。

该出版物的目的是使用众所周知的使用Python编程语言优化电驱动器和飞机的控制系统的示例,考虑在状态空间中分析和综合线性控制系统的方法。

动态优化方法的数学依据


为了进行优化,使用了振幅(MO),对称(CO)和折衷最优(KO)的标准。

解决状态空间中的优化问题时,线性平稳系统由向量矩阵方程式给出:

 dotx= fracdxdt=A cdotx+B cdotu;y=C cdotx (1)

最小控制能耗和最大速度的集成标准由功能设置:

Jx= int infty0xTQx+uTRu+2xTNudt rightarrowmin (2)

Ju= int infty0yTQy+uTRu+2yTNudt rightarrow (3)

控制律u形式为对状态变量x或对输出变量y的线性反馈:

u= pmKxu= pmKy

这样的控制使二次质量标准(2),(3)最小化。 在关系式(1)÷(3)中,Q和R分别是尺寸为[n×n]和[m×m]的对称正定矩阵; K是尺寸为[m×n]的常数系数的矩阵,其值不受限制。 如果省略输入参数N,则将其视为零。

这个问题的解决方案在状态空间中称为线性积分二次优化(LQR-optimization)问题,由以下表达式确定:

u=R1BTPx

其中矩阵P必须满足Ricatti方程:

ATP+PA+PBR1BTP+Q=0

标准(2),(3)也是矛盾的,因为第一项的实现需要一个无限大功率的源,而对于第二项,则需要无限低功率的源。 可以用以下表达式解释:

Jx= int infty0xTQxdt

这是常态 Modx 向量x,即 调节过程中其振荡的量度,因此,在瞬态中以较小的振荡获得较小的值,表达式为:

Ju= int infty0uTRudt

是用于控制的能量的量度;是对系统能量成本的惩罚。

权重矩阵Q,R和N取决于相应坐标的约束。 如果这些矩阵的任何元素等于零,则对应的坐标没有限制。 在实践中,矩阵Q,R,N的值的选择是通过专家估计,试验,误差的方法进行的,并且取决于设计工程师的经验和知识。

为了解决这些问题,使用了以下MATLAB运算符:
\开bmatrixRSE\结bmatrix=lqrABQRN\开bmatrixRSE\结bmatrix=lqryPsQRN 通过状态向量x或输出向量y最小化函数(2),(3)。

管理对象模型 Ps=ssABCD 计算结果是关于状态变量x的最佳反馈矩阵K,Ricatti方程S的解以及闭环控制系统的特征值E = eiq(A-BK)。

功能组件:

Jx=x0TP1x0Ju=x0TP2x0

其中x0是初始条件的向量; P1P2 -作为Lyapunov矩阵方程的解的未知矩阵。 使用函数可以找到它们。 P1=lyapNNQP2=lyapNNKTRKNN=A+BKT

使用Python实现动态优化方法的功能


尽管Python控制系统库[4]具有以下功能:control.lqr,control.lyap,但只有在安装了slycot模块的情况下,才能使用control.lqr,这是一个问题[5]。 在任务上下文中使用lyap函数时,优化会导致control.exception.ControlArgument:Q必须是对称矩阵错误[6]。

因此,为了实现lqr()和lyap()函数,我使用了scipy.linalg:

from numpy import* from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E P1=solve_lyapunov(NN,Ct*Q*C) 

顺便说一下,您不应该完全放弃控制,因为函数step(),pole(),ss(),tf(),feedback(),acker(),place()和其他函数都可以正常工作。

电驱动中LQR优化的示例

(示例摘自出版物[7])

考虑一个线性二次控制器的合成,该控制器满足通过矩阵在状态空间中定义的控制对象的标准(2):

A=\开bmatrix10000143.67816.667195.40201.0460 endbmatrix;B=\开bmatrix230000\结bmatrix;C=\开bmatrix100010001 endbmatrix;D=0



以下内容被视为状态变量: x1 -转换器电压,V; x2 -电动机电流,A; x3 -角速度 c1 。 这是具有HB的TP-DPT系统:发动机R nom = 30 kW; Unom = 220 V; Inom = 147 A ;; ω 0 = 169 c1 ; ω 最大值= 187 c1 ; 标称电阻力矩Mnom = 150 N * m; 起动电流的倍数= 2; 晶闸管转换器:Unom = 230 V; Uy = 10 B; Inom = 300安; 短时过电流倍数= 1.2。

解决问题时,我们接受矩阵Q对角线。 作为建模的结果,发现矩阵元素的最小值R = 84,并且 Q=[[0.01,0,0][00.01,0][0,0,0.01]] 在这种情况下,观察到发动机角速度的单调过渡过程(图1)。 在R = 840时 Q=[[0.01,0,0][00.01,0][0,0,0.01]] 瞬态过程(图2)符合MO的标准。 矩阵P1和P2的计算在x0 = [220 147 162]下进行。

列出程序(图1)。
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E #   A=matrix([[-100,0,0],[143.678, -16.667,-195.402],[0,1.046,0]]); B=matrix([[2300], [0], [0]]); C=eye(3) D=zeros([3,1]) #    R=[[84]] Q=matrix([[0.01, 0, 0],[0, 0.01, 0],[0, 0, 0.01]]); #  LQR-  [K,S,E]=lqr(A,B,Q,R) #    N=AB*K; NN=NT; w=ss(N,B,C,D) x0=matrix([[220], [147], [162]]); Ct=CT; P1=solve_lyapunov(NN,Ct*Q*C) Kt=KT; x0t=x0.T; Jx=abs((x0t*P1*x0)[0,0]) P2=solve_lyapunov(NN,Kt*R*K); Ju=abs((x0t*P2*x0)[0,0]) fig = plt.figure() plt.suptitle('      R = %s,\n\ $J_{x}$= %s; $J_{u}$ = %s,Q=[[0.01,0,0],[0, 0.01,0],[0,0,0.01]]'%(R[0][0],round(Jx,1),round(Ju,1))) ax1 = fig.add_subplot(311) y,x=step(w[0,0]) ax1.plot(x,y,"b") plt.ylabel('$x_{1}$,B') ax1.grid(True) ax2 = fig.add_subplot(312) y,x=step(w[1,0]) ax2.plot(x,y,"b") plt.ylabel('$x_{2}$,A') ax2.grid(True) ax3 = fig.add_subplot(313) y,x=step(w[2,0]) ax3.plot(x,y,"b") plt.ylabel('$x_{3}$,$C^{-1}$') ax3.grid(True) plt.xlabel(', .') plt.show() 



1个

列出程序(图2)。
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E #   A=matrix([[-100,0,0],[143.678, -16.667,-195.402],[0,1.046,0]]); B=matrix([[2300], [0], [0]]); C=eye(3) D=zeros([3,1]) #    R=[[840]] Q=matrix([[0.01, 0, 0],[0, 0.01, 0],[0, 0, 0.01]]); #  LQR-  [K,S,E]=lqr(A,B,Q,R) #    N=AB*K; NN=NT; w=ss(N,B,C,D) x0=matrix([[220], [147], [162]]); Ct=CT; P1=solve_lyapunov(NN,Ct*Q*C) Kt=KT; x0t=x0.T; Jx=abs((x0t*P1*x0)[0,0]) P2=solve_lyapunov(NN,Kt*R*K); Ju=abs((x0t*P2*x0)[0,0]) fig = plt.figure() plt.suptitle('      R = %s,\n\ $J_{x}$= %s; $J_{u}$ = %s,Q=[[0.01,0,0],[0, 0.01,0],[0,0,0.01]]'%(R[0][0],round(Jx,1),round(Ju,1))) ax1 = fig.add_subplot(311) y,x=step(w[0,0]) ax1.plot(x,y,"b") plt.ylabel('$x_{1}$,B') ax1.grid(True) ax2 = fig.add_subplot(312) y,x=step(w[1,0]) ax2.plot(x,y,"b") plt.ylabel('$x_{2}$,A') ax2.grid(True) ax3 = fig.add_subplot(313) y,x=step(w[2,0]) ax3.plot(x,y,"b") plt.ylabel('$x_{3}$,$C^{-1}$') ax3.grid(True) plt.xlabel(', .') plt.show() 



2
通过适当选择矩阵R和Q,可以将电动机的启动电流减小到等于(2-2.5)In的可接受值(图3)。 例如,R = 840且
Q =([[[0.01,0,0],[0,0.88,0],[0,0,0.01]],其值为292 A,在这些条件下的过渡过程时间为1.57 s。

列出程序(图3)。
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E #   A=matrix([[-100,0,0],[143.678, -16.667,-195.402],[0,1.046,0]]); B=matrix([[2300], [0], [0]]); C=eye(3) D=zeros([3,1]) #    R=[[840]] Q=matrix([[0.01,0,0],[0,0.88,0],[0,0,0.01]]); #  LQR-  [K,S,E]=lqr(A,B,Q,R) #    N=AB*K; NN=NT; w=ss(N,B,C,D) x0=matrix([[220], [147], [162]]); Ct=CT; P1=solve_lyapunov(NN,Ct*Q*C) Kt=KT; x0t=x0.T; Jx=abs((x0t*P1*x0)[0,0]) P2=solve_lyapunov(NN,Kt*R*K); Ju=abs((x0t*P2*x0)[0,0]) fig = plt.figure() plt.suptitle('      R = %s,\n\ $J_{x}$= %s; $J_{u}$ = %s,Q=[[0.01,0,0],[0,0.88,0],[0,0,0.01]]'%(R[0][0],round(Jx,1),round(Ju,1))) ax1 = fig.add_subplot(311) y,x=step(w[0,0]) ax1.plot(x,y,"b") plt.ylabel('$x_{1}$,B') ax1.grid(True) ax2 = fig.add_subplot(312) y,x=step(w[1,0]) ax2.plot(x,y,"b") plt.ylabel('$x_{2}$,A') ax2.grid(True) ax3 = fig.add_subplot(313) y,x=step(w[2,0]) ax3.plot(x,y,"b") plt.ylabel('$x_{3}$,$C^{-1}$') ax3.grid(True) plt.xlabel(', .') plt.show() 



图3

在所有考虑的情况下,电压和电流的反馈均为负,速度的反馈为正,根据稳定性的要求,这是不希望的。 另外,综合系统对于任务是静态的,对于负载是静态的。 因此,我们考虑在状态空间中使用附加状态变量对PI控制器进行综合 x4 -积分器的传递系数。

我们以矩阵形式表示初始信息:

A=\开bmatrix100000143.67816.667195.402001.046000010 endbmatrix;B=\开bmatrix2300000\结bmatrixC=4D=0

在R = 100,q11 = q22 = q33 = 0.001和q44 = 200的情况下,获得了与MO标准相对应的任务瞬态。图4显示了状态变量瞬态,它通过任务和负载确认了系统的静态性。

程序清单(图4)。
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt from scipy.linalg import* def lqr(A,B,Q,R): S =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*S)) E= eig(AB*K)[0] return K, S, E #   A=matrix([[-100,0,0,0],[143.678, -16.667,-195.402,0],[0,1.046,0,0] ,[0,0,1,0]]); B=matrix([[2300], [0], [0],[0]]); C=matrix([[1, 0, 0,0],[0, 1, 0,0],[0, 0, 1,0],[0, 0, 0,1]]); D=matrix([[0],[0],[0],[0]]); #    R=[[100]]; Q=matrix([[0.001, 0, 0,0],[0, 0.001, 0,0],[0, 0, 0.01,0],[0, 0,0,200]]); #  LQR-  (K,S,E)=lqr(A,B,Q,R) #    N=AB*K; NN=NT; w=ss(N,B,C,D) x0=matrix([[220], [147], [162],[0]]); Ct=CT; P1=solve_lyapunov(NN,Ct*Q*C) Kt=KT; x0t=x0.T; Jx=abs((x0t*P1*x0)[0,0]) P2=solve_lyapunov(NN,Kt*R*K); Ju=abs((x0t*P2*x0)[0,0]) fig = plt.figure(num=None, figsize=(9, 7), dpi=120, facecolor='w', edgecolor='k') plt.suptitle('      LQR -',size=14) ax1 = fig.add_subplot(411) y,x=step(w[0,0]) ax1.plot(x,y,"b") plt.ylabel('out(1)') ax1.grid(True) ax2 = fig.add_subplot(412) y,x=step(w[1,0]) ax2.plot(x,y,"b") plt.ylabel('out(2)') ax2.grid(True) ax3 = fig.add_subplot(413) y,x=step(w[2,0]) ax3.plot(x,y,"b") plt.ylabel('out(3)') ax3.grid(True) ax4 = fig.add_subplot(414) y,x=step(w[3,0]) ax4.plot(x,y,"b") plt.ylabel('out(4)') ax4.grid(True) plt.xlabel(', .') plt.show() 



4
为了确定矩阵K,Python控制系统库具有两个函数K = acker(A,B,s)和K = place(A,B,s),其中s是向量-闭环控制系统传递函数的期望极点所在的行。 第一条命令仅可用于u中n≤5且u中有一个输入的系统。 第二个没有这种限制,但是极点的多重性不应该超过矩阵B的秩。清单和(图5)中显示了使用acker()运算符的示例。

程序清单(图5)。
 # -*- coding: utf8 -*- from numpy import* from control.matlab import * from matplotlib import pyplot as plt i=(-1)**0.5 A=matrix([[-100,0,0],[143.678, -16.667,-195.402],[0,1.046,0]]); B=matrix([[2300], [0], [0]]); C=eye(3) D=zeros([3,1]) p=[[-9.71+14.97*i -9.71-14.97*i -15.39 -99.72]]; k=acker(A,B,p) H=AB*k; Wss=ss(H,B,C,D); step(Wss) w=tf(Wss) fig = plt.figure() plt.suptitle('   acker()') ax1 = fig.add_subplot(311) y,x=step(w[0,0]) ax1.plot(x,y,"b") ax1.grid(True) ax2 = fig.add_subplot(312) y,x=step(w[1,0]) ax2.plot(x,y,"b") ax2.grid(True) ax3 = fig.add_subplot(313) y,x=step(w[2,0]) ax3.plot(x,y,"b") ax3.grid(True) plt.xlabel(', .') plt.show() 




5

结论


考虑到新lqr()和lyap()函数的使用,该出版物中提出的LQR驱动优化的实现将使我们在研究控制理论的相应部分时放弃使用许可的MATLAB程序。

飞机(LA)中LQR优化的示例

(示例取自Astrom和Mruray的出版物,第5章[8],并由文章的作者最终完成,以取代lqr()函数并将术语引入公认的术语)

理论部分,简要的理论,LQR优化已经在上面考虑过了,所以让我们继续进行代码分析并给出相应的注释:

所需的库和LQR控制器功能。

 from scipy.linalg import* # scipy   lqr from numpy import * # numpy    from matplotlib.pyplot import * #     from control.matlab import * #  control..ss()  control..step() #     lqr() m = 4; #    . J = 0.0475; #        ***2 #́-       #  . r = 0.25; #     . g = 9.8; #     /**2 c = 0.05; #   () def lqr(A,B,Q,R): X =matrix(solve_continuous_are(A, B, Q, R)) K = matrix(inv(R)*(BT*X)) E= eig(AB*K)[0] return X, K, E 

初始条件和基本矩阵A,B,C,D模型。

 xe = [0, 0, 0, 0, 0, 0];#   x   () ue = [0, m*g]; #   #   A = matrix( [[ 0, 0, 0, 1, 0, 0], [ 0, 0, 0, 0, 1, 0], [ 0, 0, 0, 0, 0, 1], [ 0, 0, (-ue[0]*sin(xe[2]) - ue[1]*cos(xe[2]))/m, -c/m, 0, 0], [ 0, 0, (ue[0]*cos(xe[2]) - ue[1]*sin(xe[2]))/m, 0, -c/m, 0], [ 0, 0, 0, 0, 0, 0 ]]) #   B = matrix( [[0, 0], [0, 0], [0, 0], [cos(xe[2])/m, -sin(xe[2])/m], [sin(xe[2])/m, cos(xe[2])/m], [r/J, 0]]) #   C = matrix([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]]) D = matrix([[0, 0], [0, 0]]) 

我们建立与xy位置的逐步变化相对应的输入和输出。 向量xd和yd对应于系统中所需的稳态。 矩阵Cx和Cy是模型的相应输出。 为了描述系统的动力学,我们使用矢量矩阵方程组:

\左\ {\开始{矩阵} xdot = Ax + Bu => xdot =(A-BK)x + xd \\ u = -K(x-xd)\\ y = Cx \\ \结束{matrix} \对。



对于K \ cdot xd和K \ cdot xd作为输入向量,可以使用step()函数对闭环的动力学进行建模,其中:

 xd = matrix([[1], [0], [0], [0], [0], [0]]); yd = matrix([[0], [1], [0], [0], [0], [0]]); 

当前控件0.8.1库仅支持用于创建反馈控制器的SISO Matlab函数的一部分,因此,要从控制器读取数据,因此有必要为横向-x和垂直-y动态创建索引向量lat(),alt()。

 lat = (0,2,3,5); alt = (1,4); #    Ax = (A[lat, :])[:, lat]; Bx = B[lat, 0]; Cx = C[0, lat]; Dx = D[0, 0]; Ay = (A[alt, :])[:, alt]; By = B[alt, 1]; Cy = C[1, alt]; Dy = D[1, 1]; #      K Qx1 = diag([1, 1, 1, 1, 1, 1]); Qu1a = diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1a) K1a = matrix(K); #      x H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx); (Yx, Tx) = step(H1ax, T=linspace(0,10,100)); #     y H1ay = ss(Ay - By*K1a[1,alt], By*K1a[1,alt]*yd[alt,:], Cy, Dy); (Yy, Ty) = step(H1ay, T=linspace(0,10,100)); 

图形按顺序显示,每个任务一个。

 figure(); title("   x  y") plot(Tx.T, Yx.T, '-', Ty.T, Yy.T, '--'); plot([0, 10], [1, 1], 'k-'); axis([0, 10, -0.1, 1.4]); ylabel(' '); xlabel(''); legend(('x', 'y'), loc='lower right'); grid() show() 

图表:



输入动作幅度的影响
 #     Qu1a = diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1a) K1a = matrix(K); H1ax = ss(Ax - Bx*K1a[0,lat], Bx*K1a[0,lat]*xd[lat,:], Cx, Dx); Qu1b = (40**2)*diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1b) K1b = matrix(K); H1bx = ss(Ax - Bx*K1b[0,lat], Bx*K1b[0,lat]*xd[lat,:],Cx, Dx); Qu1c = (200**2)*diag([1, 1]); X,K,E =lqr(A, B, Qx1, Qu1c) K1c = matrix(K); H1cx = ss(Ax - Bx*K1c[0,lat], Bx*K1c[0,lat]*xd[lat,:],Cx, Dx); figure(); [Y1, T1] = step(H1ax, T=linspace(0,10,100)); [Y2, T2] = step(H1bx, T=linspace(0,10,100)); [Y3, T3] = step(H1cx, T=linspace(0,10,100)); plot(T1.T, Y1.T, 'r-'); plot(T2.T, Y2.T, 'b-'); plot(T3.T, Y3.T, 'g-'); plot([0 ,10], [1, 1], 'k-'); axis([0, 10, -0.1, 1.4]); title("   ") legend(('1', '$1,6\cdot 10^{3}$','$4\cdot 10^{4}$'), loc='lower right'); text(5.3, 0.4, 'rho'); ylabel(' '); xlabel(''); grid() show() 


图表:



瞬态响应
 #  Qx    Qx2 = CT * C; Qu2 = 0.1 * diag([1, 1]); X,K,E =lqr(A, B, Qx2, Qu2) K2 = matrix(K); H2x = ss(Ax - Bx*K2[0,lat], Bx*K2[0,lat]*xd[lat,:], Cx, Dx); H2y = ss(Ay - By*K2[1,alt], By*K2[1,alt]*yd[alt,:], Cy, Dy); figure(); [Y2x, T2x] = step(H2x, T=linspace(0,10,100)); [Y2y, T2y] = step(H2y, T=linspace(0,10,100)); plot(T2x.T, Y2x.T, T2y.T, Y2y.T) plot([0, 10], [1, 1], 'k-'); axis([0, 10, -0.1, 1.4]); title(" ") ylabel(' '); xlabel(''); legend(('x', 'y'), loc='lower right'); grid() show() 


图表:



基于物理的瞬态响应
 #    #       #   1 .     10 . Qx3 = diag([100, 10, 2*pi/5, 0, 0, 0]); Qu3 = 0.1 * diag([1, 10]); X,K,E =lqr(A, B, Qx3, Qu3) K3 = matrix(K); H3x = ss(Ax - Bx*K3[0,lat], Bx*K3[0,lat]*xd[lat,:], Cx, Dx); H3y = ss(Ay - By*K3[1,alt], By*K3[1,alt]*yd[alt,:], Cy, Dy); figure(); [Y3x, T3x] = step(H3x, T=linspace(0,10,100)); [Y3y, T3y] = step(H3y, T=linspace(0,10,100)); plot(T3x.T, Y3x.T, T3y.T, Y3y.T) axis([0, 10, -0.1, 1.4]); title("   ") ylabel(' '); xlabel(''); legend(('x', 'y'), loc='lower right'); grid() show() 


图表:



结论:


考虑到对新lqr()函数的使用,在出版物中介绍的飞机中LQR优化的实现允许我们在研究控制理论的相应部分时放弃使用许可程序MATLAB。

参考文献


1. 手指上的数学:线性二次控制器。

2. 设计最优控制系统问题中的状态空间。

3. 使用Python控制系统库设计自动控制系统。

4. Python控制系统库。

5. Python-无法将slycot模块导入spyder(RuntimeError和ImportError)。

6. Python控制系统库。

7. 电动驱动器中最佳控制和lqr优化的标准。

Source: https://habr.com/ru/post/zh-CN441582/


All Articles