Python解决控制问题系列之二：线性连续系统最优控制问题

文章目录

Python解决控制问题系列之二：线性连续系统最优控制问题
- 1. 前言
- 2. 线性系统问题描述
- 3. Python 编程
- - 3.1 仿真状态模型构造
  - 3.2 黎卡提方程求解
  - 3.3 构造反馈控制器的I/O系统描述
  - 3.4 闭环系统构造
- 4 结语

1. 前言

线性系统是控制问题的最常见状态空间表达式模型，也是各类物理运动系统在平衡点处线性化后的标准模型。本文是在研究如何利用python编程复现Lewis¹科研团队在2009年Automatica发表的一篇有关“reinforcement learning control”的算法时，衍生出来的基础线性系统模型构造、代数黎卡提方程（ARE）求解以及状态反馈闭环系统的轨迹数值解求解实现方案。

2. 线性系统问题描述

线性系统的模型如下式：
x˙=Ax+Bu(1)\dot x = Ax + Bu\quad (1) x˙=Ax+Bu(1)
最优控制的价值函数为无穷时间尺度，要求找到最优控制u∗u^*u∗,使得下列价值函数取得最小值，同时保证系统稳定。
J=∫t0∞xT(τ)Qx(τ)+uT(τ)Ru(τ)dτ(2)J = \int^{\infin}_{t_0} x^T(\tau)Qx(\tau)+u^T(\tau)Ru(\tau)d\tau \quad (2) J=∫t0∞xT(τ)Qx(τ)+uT(τ)Ru(τ)dτ(2)
通过最优控制的基本理论分析可知，为例保证（2）取得最小值，必须保证从代数黎卡提方程中求得对称正定阵点P，即：
ATP+PA−PBR−1BTP+Q=0(3)A^TP+PA-PBR^{-1}B^TP + Q = 0 \quad (3) ATP+PA−PBR−1BTP+Q=0(3)
此时最优控制律u*可以表示为：
u∗(t)=−Kx(t)=−R−1BTPx(t)(4)u^*(t) = - K x(t) = -R^{-1}B^TPx(t)\quad (4) u∗(t)=−Kx(t)=−R−1BTPx(t)(4)
从上述（4）式可知，问题的核心就是求解代数黎卡提方程ARE和构建闭环控制系统。

3. Python 编程

由于scipy.integrate.solve_ivp()函数只适用于一维的微分方程数值求解，这样当面对多维的线性系统状态方程时，就必续将系统状态方程人为转化为每个变量的一维微分方程，非常不方便，也无法利用后续的矩阵特性求解。因此，采用python-control功能包中的若干函数API来完成状态方程的构造、控制器的构造、输入输出轨迹响应数值求解的构造。

3.1 仿真状态模型构造

仿真系统来源于文献【1】中提到的能量系统在关键操作点的简化数学模型，即：
A=(−0.066511.5000−2.52.50−9.50−13.736−13.7360.6000)B=(0013.7360)TA = \begin{pmatrix} -0.0665 & 11.5 & 0 & 0 \\ 0 & -2.5 & 2.5 & 0 \\ -9.5 & 0 & -13.736 & -13.736\\ 0.6 & 0 & 0 & 0 \end{pmatrix}\\ B = \begin{pmatrix}0 & 0 & 13.736 & 0\end{pmatrix}^T A=⎝⎜⎜⎛−0.06650−9.50.611.5−2.50002.5−13.736000−13.7360⎠⎟⎟⎞B=(0013.7360)T
通过control.ss函数实现状态空间表达式的构建，将C设置为单位矩阵，保证输出为全状态，D设置为0矩阵（不考虑输入对输出的影响）。当状态方程构建完成后，再加线性系统状态方程转化为I/O系统描述，方便后续调用control.input_output_response API接口函数实现轨迹的数值解求解。具体的代码为：

# 引入必要的功能包和API接口
import numpy as np
import control as ct
import matplotlib.pyplot as plt
# 科学计数法, 小数精度控制在小数点后3位
np.set_printoptions(precision=4,suppress=True)
#构造线性系统 \dot x = Ax+Bu
A = np.array([[-0.0665, 11.5, 0, 0],[0, -2.5, 2.5, 0],[-9.5, 0, -13.736, -13.736],[0.6, 0, 0, 0]
])
B = np.array([[0], [0], [13.736], [0]])
C = np.diag([1.0, 1.0, 1.0, 1.0])
D = np.zeros((4, 1))
#构造出状态方程
sys = ct.ss(A, B, C, D)
# 将状态方程转化为I/O系统描述，定义输入名为controller,输出为四个状态。
state_sys = ct.LinearIOSystem(sys, name='state_sys', inputs=('controller'), outputs=('x[0]', 'x[1]','x[2]','x[3]'), states=4 )

构造出的state_sys I/O系统画成方框图表示即为：

3.2 黎卡提方程求解

control功能包提供了一个函数接口control.care来方便计算标准的代数黎卡提方程，如（3）式。具体的API接口使用说明再次不在赘述，可以参考control包的使用说明。与文献【1】中的设置相同，将对称阵Q和R均设置为单位矩阵，从而通过求解代数黎卡提方程，得到三个返回值，分别是解P，闭环系统的特征值L，状态反馈增益K，这就为后续构造状态反馈控制器做好了准备。具体代码为：

# 初始化Q
Q = np.diag([1.0, 1.0, 1.0 , 1.0])
# 解标准的代数黎卡提方程，R=None时，默认R为单位阵
P, L, K = ct.care(A, B, Q, R=None)

3.3 构造反馈控制器的I/O系统描述

接下去就是将基于全状态反馈的控制器构造成I/O系统，方便和之前构造的state_sys I/O系统进行反馈链接，实现初值状态下闭环系统的构建。控制器不涉及到微分方程，所以在使用control.NonlinearIOSystem接口时，第一项参数填为None，第二项参数作为I/O系统的输出项，即反馈控制器的表达式 u = -kx。值得注意的是，在设置输入端口时，需要多设置一个控制输入（这个控制输入在实际运行中不起任何作用，设置的原因是因为在使用control.input_output_response时必须设置初始闭环系统的输入，在这里你可以把它当作一个启动触发信号，对状态的轨迹无任何影响）。具体代码如下：

# 定义反馈控制器ctrl = -Kx，其中参数u为5维列向量，前4维为四个状态，最后1维为启动触发信号ref
def control_output(t, x, u, params={'control gain': np.zeros((1, 4))}):k = params['control gain']  #默认K为0,开环控制x = u[:4]ctrl = -k @ x return ctrl
# 构造控制器 I/O 系统
control_sys = ct.NonlinearIOSystem(None, control_output, params={'control gain': K}, name='control_sys', inputs=('x[0]', 'x[1]','x[2]','x[3]', 'ref'),outputs = 'u')

control_sys 构造成的I/O系统化程方框图为：

3.4 闭环系统构造

通过control.interconnect API接口实现将控制对象state_sys I/O系统的输入输出端与控制器control_sys I/O 系统的输入输出端相连接，构成闭环控制系统。即：

对于闭环系统closed_loop系统来说，它的输入即ref信号，作为仿真启动信号，取值为0即可，对具体轨迹变化无影响。输出为所需的四位状态轨迹x[0]，…, x[3]。系统状态的初始值取为：
x0=[0,0.1,0,0]x_0 = [0, 0.1, 0, 0] x0=[0,0.1,0,0]
具体的代码为：

#构造闭环系统
closed_loop = ct.interconnect((state_sys, control_sys), connections=(['state_sys.controller', 'control_sys.u'],['control_sys.x[0]', 'state_sys.x[0]'],['control_sys.x[1]', 'state_sys.x[1]'],['control_sys.x[2]', 'state_sys.x[2]'],['control_sys.x[3]', 'state_sys.x[3]']),inplist=['control_sys.ref'],inputs=['start'], outlist=['state_sys.x[0]','state_sys.x[1]','state_sys.x[2]','state_sys.x[3]'],outputs=['x[0]','x[1]','x[2]','x[3]']
)#计算轨迹响应
t = np.array([i * 0.5 for i in range(100)]) #每隔0.5取一个点进行数值计算
tp = np.linspace(0, 50, 200)
tout, yout = ct.input_output_response(sys, tp, 0, [0, 0.1, 0, 0], t_eval = t) #开环控制的轨迹,按照t的采用分布取值
t1, y1 = ct.input_output_response(closed_loop, tp, 0, [0, 0.1, 0, 0], t_eval = t) #闭环控制的轨迹
plt.plot(tout, yout[1,:])
plt.plot(t1, y1[1, :])
print(t1)

最终得到的效果图为：

蓝色线为x[1]状态的开环控制轨迹变化情况，橙色线为x[1]状态在闭环最优控制下轨迹变化情况，收敛到平衡点的速度明显加快。

4 结语

查了很多资料，网上显有python实现线性系统控制的程序范例，大都都是matlab实现。后续，如果攻破了文献【1】的算法复现的话，还会共享出来。

D. Vrabie, O. Pastravanu, M. Abu-khalaf, F.L. Lewis. Adaptive optimal control for continuous-time linear systems based on policy iteration. Automatica, 2009(45): 477-484 ↩︎