我正在尝试解决受动态(h(x, x', u) = 0)约束的非线性最优控制问题。
给定:
f(x) = (u(t) - u(0)(t))^2 # u0(t) 是提供给系统的初始输入
h(x) = y'(t) - 积分(sqrt (u(t))*y(t) + y(t)) = 0 # 一个非线性微分方程
-2 < y(t) < 10 # 系统状态限制在这个范围内
-2 < u(t) < 10 # 系统状态限制在这个范围
u0(t) # 将被定义为任意分段线性函数
我尝试使用 openopt 和 scipy 将问题转换为 python 代码:
import numpy as np
from scipy.integrate import *
from openopt import NLP
import matplotlib.pyplot as plt
from operator import and_
N = 15*4
y0 = 10
t0 = 0
tf = 10
lb, ub = np.ones(2)*-2, np.ones(2)*10
t = np.linspace(t0, tf, N)
u0 = np.piecewise(t, [t < 3, and_(3 <= t, t < 6), 6 <= t], [2, lambda t: t - 3, lambda t: -t + 9])
p = np.empty(N, dtype=np.object)
r = np.empty(N, dtype=np.object)
y = np.empty(N, dtype=np.object)
u = np.empty(N, dtype=np.object)
ff = np.empty(N, dtype=np.object)
for i in range(N):
t = np.linspace(t0, tf, N)
b, a = t[i], t[i - 1]
integrand = lambda t, u1, y1 : np.sqrt(u1)*y1 + y1
integral = lambda u1, y1 : fixed_quad(integrand, a, b, args=(u1, y1))[0]
f = lambda x1: ((x1[1] - u0[i])**2).sum()
h = lambda x1: x1[0] - y0 - integral(x1[0], x1[1])
p[i] = NLP(f, (y0, u0[i]), h=h, lb=lb, ub=ub)
r[i] = p[i].solve('scipy_slsqp')
y0 = r[i].xf[0]
y[i] = r[i].xf[0]
u[i] = r[i].xf[1]
ff[i] = r[i].ff
figure1 = plt.figure()
axis1 = figure1.add_subplot(311)
plt.plot(u0)
axis2 = figure1.add_subplot(312)
plt.plot(u)
axis2 = figure1.add_subplot(313)
plt.plot(y)
plt.show()
现在的问题是,使用像y0 = 10这样的正初始 y0 运行代码,代码将产生令人满意的结果。但是给y0 = 0或负一个y0 = -1,nlp 问题将是有缺陷的,说:
“没有获得可行的解决方案(1 个约束等于 NaN,MaxResidual = 0,objFunc = nan) ”
另外,考虑到分段线性初始 u0,如果在 函数的第一个范围内放置 0 以外的任何数字t < 3
,则表示:
u0 = np.piecewise(t, [t < 3, and_(3 <= t, t < 6), 6 <= t], [2, lambda t: t - 3, lambda t: -t + 9])
而不是:
u0 = np.piecewise(t, [t < 3, and_(3 <= t, t < 6), 6 <= t], [0, lambda t: t - 3, lambda t: -t + 9])
这将再次导致相同的错误。
有任何想法吗 ?
提前致谢。