打靶法（Shooting Method）在Python中的实现（SciPy库的使用）

本文为笔者在学习最优控制问题的传统数值解法——“打靶法”——过程中所做的尝试，因为对Matlab语法不熟悉所以想尝试在Python中实现这些操作，参考了其他大佬们的代码学习到了很多，也希望能给其他人带来帮助。将不定时更新。因为是编程小白所以代码可能会出现错误或者是写得非常不好看，希望大家及时指出交流。所有代码均在Jupyter上编译运行，请大家根据自己搭建的环境更改。

单步直接打靶法Direct Single Step Shooting Method (DSSSM)

参考的最爱大盘鸡大佬的代码：[https://blog.csdn.net/Ruins_LEE/article/details/125680062?spm=1001.2014.3001.5502]

该方法的要点是，离散控制变量，将目标函数、状态变量、约束条件、边界条件等转化为控制变量的函数，进行数值求解。

DSSSM的Python实现

算例

该算例选自《最优化与最优控制》第2版第257页例13.1。

设由状态方程及初始条件$\dot{x}=-x^2+u$，$x_0=10$，性能指标$J(u)=0.5{\int }_0^1(x^2+u^2)dt$，，求解最优控制使$J$为极小。

import numpy as np
from scipy.optimize import minimize
from scipy.interpolate import interp1d
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

"""Parameters Setting"""
# 初始参数设置  【使用词典更方便传递多个参数】
system = {'t0': 0, 'tf': 1, 'x_0': 10}
# t0:起始时间  tf:终止时间  x_0:状态变量初值
system['N'] = 20 # 打靶点参数设置
system['t'] = np.linspace(system['t0'], system['tf'], system['N']) #离散化的时间

# 离散控制变量
#system['u_index'] = np.arange(1, system['N'] + 1)

# 先猜测一组控制变量
u_0 = - 0.5 * np.ones(system['N'])
system['u'] = u_0

"""Subfunctions Setting"""
# 定义目标函数
def objective_J(u, p):
    # 使用离散的控制变量计算离散的状态变量
    #使用ode解算器 计算状态变量x
    sol = solve_ivp(lambda t, x: derivative_dx(t, x, p), [p['t0'], p['tf']], [p['x_0']]
                    , t_eval = p['t'])
    x = sol.y.flatten()
    # 计算目标函数值
    L = x ** 2 + u ** 2
    #梯形规则进行积分
    int_L = 0.5 * np.trapz(L, p['t']) 
    return int_L
# 定义状态方程
def derivative_dx(t, x, p):
    # interp1d进行插值
    u_interp = interp1d(p['t'], p['u'], fill_value = "extrapolate")
    u = u_interp(t)
    dxdt = - x ** 2 + u
    return dxdt

"""Solving Algorithm"""
# 求解控制变量
result = minimize(objective_J, u_0, args=(system,), method='SLSQP', options={'disp': True})
# 更新词典里的控制变量
system['u'] = result.x

# 使用ode解算器计算控制变量对应的状态变量
sol = solve_ivp(lambda t, x: derivative_dx(t, x, system), [system['t0'], system['tf']], [system['x_0']]
                , t_eval=system['t'])
system['x'] = sol.y.flatten()

plt.figure(figsize=(10, 5))
plt.plot(system['t'], system['x'], 'x-', label='State $x(t)$', linewidth=1.5)
plt.plot(system['t'], system['u'], 'x-', label='Control $u(t)$', linewidth=1.5)
plt.xlabel('Time')
plt.ylabel('State & Control')
plt.title('State and Control versus Time')
plt.legend()
plt.grid(True)

plt.savefig('DSSSM.png')
plt.show()

运行结果