python 求解哈密顿最优控制程序例子

在Python中，我们可以使用一些优化库如`scipy.optimize`或更专业的数学优化库如`pyomo`、`cvxopt`来解决哈密顿雅可比贝尔曼（HJB）方程或者动态规划问题，以求解最优控制问题。这里举一个简单的例子，我们假设有一个经典的连续时间线性二次型最优化问题：

import numpy as np
from scipy.optimize import minimize

# 定义状态空间矩阵 A 和输入矩阵 B
A = np.array([[0, 1], [-1, -2]])
B = np.array([[0], [1]])

# 目标函数Q(t,x,u)
def qfun(x, u):
    return x[0]**2 + x[1]**2 + u**2

# 边界条件和目标函数
def cost_function(x, u, T):
    return integral(qfun, (0, T), args=(x, u))

def integral(f, a, b, *args):
    # 使用辛普森法则计算积分
    dt = 0.01
    integral_sum = f(a, *args) + f(b, *args)
    for i in range(1, int((b-a)/dt)):
        integral_sum += 4 * f(a + i*dt, *args) + 2 * f(a + (i+1)*dt, *args)
    return dt * integral_sum / 3

# 设置初始状态、最终时间和其他约束
x0 = np.array([1, 0])
T = 5
u_ub = np.inf
u_lb = -np.inf

# 构建优化问题
cons = ({'type': 'ineq', 'fun': lambda x, _: u_lb <= x[1]},
        {'type': 'ineq', 'fun': lambda x, _: x[1] <= u_ub})
bnds = [(None, None)] * len(x0)

# 求解最优化问题
res = minimize(cost_function, x0, method='SLSQP', bounds=bnds, constraints=cons,
               options={'disp': True})

# 输出结果
print("Optimal state:", res.x)
print("Optimal control:", res.fun)