python解决控制问题系例之一：解决微分方程-状态方程求解作图问题

python解决控制问题系例之一：解决微分方程-状态方程求解作图问题

1. 利用scipy.integrate功能包的solve_ivp求解方案

以求解lotka-Volterra 系统为例来演示如何使用slove_ivp。lotka-Volterra系统的动态模型可以描述为：
$$\begin{gathered} \dot{x}=ax-bxy, \\ \dot{y}=-cy+dxy \end{gathered}$$
在使用solve_ivp前，需要先建立系统模型函数lotkavolterra，函数的输入为自变量时间t、输入状态向量z=(x, y)以及系统的参数a,b,c,d。函数的返回值为一个列表，列表中每一项为一行微分方程的更新值。该部分具体的代码为：

# 引入必要的功能包
from scipy.integrate import solve_ivp
import numpy as np
import matplot.pyplot as plt

def lotkavolterra(t, z, a, b, c, d):
    x, y = z
    return [a * x - b * x * y, -c * y + d * x * y]

调用solve_ivp来生成微分方程的数值解，包含的参数包括：系统模型函数lotkavolterra, 时间t的范围，状态初值，以及参数a、b、c、d的具体值。具体代码为：

sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1), dense_output=True)
t = np.linspace(0, 15, 300)
z = sol.sol(t)

最后使用matplotlib功能包中的plot功能，画出系统状态关于时间t的轨迹。

plt.plot(t, z.T)
plt.legend(['x', 'y'], loc='best', shadow=True)
plt.xlabel('t')
plt.title('Lotka-Volterra System')
plt.show()

最终效果图：

2. 利用python-control功能包的input_output_response求解方案

此方案中需要用到control功能包。control功能包是包含了很多控制系统中常用的建模、分析和类matlab函数使用方法，对于求解微分方程构成的系统轨迹非常方便。

以功能包中的案例 predator-prey为例，讲解微分方程数值解求解方法。首先引入必要的功能包，包括numpy、control和matplotlib.pyplot, 代码为：

import control as ct
import numpy as np
import matplotlib.pyplot as plt

构建predator_prey的模型，其微分方程为：
$$\begin{gathered} \dot{H}=\frac{(r+u)H}{1-H/K}-\frac{aHL}{c+H}, \\ \dot{L}=b\frac{aHL}{c+H}- \end{gathered}$$