Python科学计算包在最优控制问题中的应用

Python科学计算包在最优控制问题中的应用

岳玉环

嘉兴职业技术学院 智能制造学院 浙江省嘉兴市 314036

摘　要：在控制工程中，常采用微分方程（组）来描述系统，但利用解析方法求解微分方程（组）的计算难度较大，更多采用数值方法进行求解。以非线性弹簧-质量系统的最优控制问题为例，展示了利用开源Python 的科学计算包Scipy可以方便、直观、高效的实现非线性微分方程的数值求解，对Python在自动控制领域的推广应用有很好的借鉴意义。

关键词：Python；科学计算包Scipy；微分方程（组）；数值解；最优控制

中图分类号：TP391

投稿投稿

0引言

众所周知，在自动控制工程中，经常采用微分方程（组）来描述系统特征。然而用解析方法来求解微分方程（组）是非常困难的，为此在实际的控制工程中，更多的采用数值求解方法[1]。能够实现数值运算有多种专用软件，如Mathematica, Matlab, Python等。而Python语言因其开源的特性及简洁、易读和可扩展性，从众多的计算软件中脱颖而出，已成为最受欢迎的程序设计语言。

Python 专用的科学计算包Scipy 可以方便的实现各类实际工程涉及到的数值计算功能[2]。在Scipy 科学计算包中包含诸如Scipy.integrate（积分程序）、Scipy.optimize（优化）、Scipy.interpolate （插值）、Scipy.fftpack （傅里叶变换）、Scipy.signal （信号处理）、Scipy.linalg （线性代数程序）等十几个功能子包。每个子包中提供了若干函数和类用于解决各自对应领域的数值计算问题。本文以Scipy.integrate为例，介绍其中的solve_bvp函数在实际控制工程中的应用方法，以此展示Python科学计算包Scipy简洁、直观、高效的实现各类型微分方程（组）的数值求解，对Python在自动控制领域的推广应用有很好的借鉴意义。

1非线性弹簧-质量系统的最优控制问题

本节中通过非线性弹簧-质量系统的最优控制问题展示solve_bvp函数求解微分方程（组）的边值问题的一般过程。

1.1 非线性弹簧-质量系统方程及问题描述

对于实际的弹簧-质量系统更多的情况下是非线性的。本文选择一个非线性弹簧-质量系统[3]，其方程表达式为

⑴

其中为物体偏离原点的位移，此处的原点定义为，也就是弹簧没有受力时的自由状态物体所在的位置。为与弹簧连接的物体质量。和分别为线性和非线性弹簧系数。 为该系统的控制输入信号。我们期望通过选择合适的控制信号使得该非线性弹簧-质量系统从已知的初始状态能够以某种最优方式（如能量消耗最少）到达原点处()并稳定下来。这类最优控制问题可以借助于Pontryagin极大原理转换为两点边值问题（bvp）进行求解。

首先，将原系统方程(1)转化为一阶微分方程的形式。为了简单起见，我们选择，则可以得到

假设,得到一阶微分方程组

进一步写作反馈控制的形式，其中，此处以能量方式构造消耗函数（Cost Function）如下：

根据Pontryagin极大原理[4]，从初始状态出发的局部最优控制的轨道满足如下方程：

⑵

并满足边值条件，其中的

为对偶曲线。求解一阶微分方程组(2})得到和。利用即可以得到对应在初始状态时，使得性能指标J最小的最优控制信号[5]。

1.2 利用Scipy.integrate.solve_bvp求解最优控制信号

该函数形式为：scipy.integrate.solve_bvp(fun,bc,x,y,p=None,S=None,fun_jac=None,bc_jac=None,tol=0.001,max_nodes=1000,verbose=0,bc_tol=None)，其中各参数说明详见网站资料[6]。下面是应用Scipy.integrate.solve_bvp函数进行最优控制信号求解的主要过程。

1.2.1 参数定义

（1）将要求解的方程(2)定义为方向场函数fun(x,y)。具体程序语句：

def fun(x,y):

dot_y_0 = y[1]

dot_y_1 = - y[0] - xi*y[0]**3-y[3]

dot_y_2 = y[3] - y[0] + 3*xi*y[0]**2*y[3]

dot_y_3 = - y[2]-y[1]

dot_y_0_1 = np.vstack((dot_y_0,dot_y_1))

dot_y_2_3 = np.vstack((dot_y_2,dot_y_3))

return np.vstack((dot_y_0_1,dot_y_2_3))

（2）给出边界条件。在本例中假设系统的初始状态为，最终到达原点并稳定下来，即另一边界条件为，将该边界条件转化为。在进行数值计算时常将条件具体化为某个较大的数值，如本例中如采用20来替换。具体程序语句为：

def bc(ya,yb):

return np.vstack(([ya[0:2] - np.array([3.0,-0.5]),yb[2:4]]))

（3）构造自变量的网格点，点数是任意的，但区间应该与前面的边界条件一致的，即，如前所述本例中自变量边界为到。具体程序语句为：

x = np.linspace(0, 20, 5)

（4）给出一组在自变量网格点处函数的初始猜测值，该值是任意的。本例中给出的初始猜测值为。具体程序语句为：

y_cc = np.zeros((4,x.size))

（5）其他参数采用默认数值即可。

1.2.2 求解最优控制信号

（1）直接调用Scipy.integrate.solve_bvp即可求解方程(2)。在函数求解后返回的数组.sol就是对应的和。具体程序语句为：

res_cc = solve_bvp(fun, bc , x , y_cc)

（2）求解得到最优控制信号，具体程序语句为：

x_plot = np.linspace(0,20,100)

p_plot_cc = res_cc.sol(x_plot)[2:4,]

u_pot_cc = - np.dot(g, p_plot_cc)

（3）对运算结果进行可视化处理[7]，绘制位移曲线和最优控制信号曲线。为了对比控制效果，分别给出了两种初始状态下的系统响应及最优控制信号。图1为对应初始状态的位移曲线和最优控制信号曲线，图2为对应初始状态的位移曲线和最优控制信号曲线。

图1：对应初始状态曲线

图2：对应初始状态曲线

2结束语

通过在非线性弹簧-质量系统的最优控制问题中的具体应用，展示了scipy科学计算包中关于边值问题进行数值计算的强大功能。除如前所述，scipy科学计算包中还提供其他满足各类工程中的实际应用问题的数值计算工具，为解决自动控制工程领域中遇到的数值求解问题提供了更加低成本，高效率的新思路和新方法，对Python在自动控制领域的推广应用有很好的借鉴意义。

3参考文献

[1]K. J. Astrom and R. M. Murray, Feedback Systems: An Introduction For Scientists And Engineers, Princeton University Press, Princeton, NJ, USA, 2008.

[2]Python科学计算[M]. 清华大学出版社 , 张若愚, 2016.

[3]Sakamoto, N, van, et al. Analytical Approximation Methods for the Stabilizing Solution of the Hamilton–Jacobi Equation[J]. IEEE Transactions on Automatic Control, 2008, 53(10):2335-2350.

[4]Hassan K Khalil. Nonlinear systems [M]. Upper Saddle River, 2002.

[5]Daniel Liberzon. Calculus of variations and optimal control theory, Princeton University Press, 2012.

[6]https://docs.scipy.org/doc/scipy/reference/index.html

[7]Jesse M.Kinder. Python物理建模初学者指南[M]. 盖磊 译.人民邮电出版社,2017.