我试图解决一个受动态约束的非线性最优控制问题(h(x,x ',u)= 0)。
给定:
f(x)=(u(t)- u(0)(t))^2#u 0(t)是提供给系统的初始输入
*h(x)= y '(t)-积分(sqrt(u(t))y(t)+ y(t))= 0#一个非线性微分方程
-2〈y(t)〈10#系统状态有界于此范围
-2〈u(t)〈10#系统状态有界于此范围
u 0(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()
现在的问题是,运行一个正的初始y 0的代码,如y 0 = 10,代码将得到满意的结果。但给定y 0 = 0或一个负的y 0 = -1,nlp问题将是有缺陷的,即:
“未获得可行解(1个约束等于NaN,MaxResidual = 0,objFunc = nan)”
同样,考虑分段线性初始值u 0,* 如果 * 在t < 3
处的函数的第一个值域中放置除0以外的任何数字,则意味着: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])
这将再次导致相同的错误。
有什么想法吗?
先谢谢你。
2条答案
按热度按时间jjjwad0x1#
我的第一个想法是,你似乎在解决一个二维的最优控制问题,就好像它是一个一维的问题。
约束动力学$h(x,x ',t)$实际上是一个二阶常微分方程。
(t)- sqrt(u(t))*y(t)+ y(t))= 0
从这个开始,我将以标准的方式将我的系统重新定义为二维一阶系统。
我的第二个想法是,你似乎是在每个时间步对$u(t)$进行独立优化,而问题是对$u(.)$(整个函数)进行全局优化。因此,如果有什么不同的话,对NLP的调用应该在for循环之外......
有专门的Optimal Control开放源代码工具箱:
Python中有JModellica:http://www.jmodelica.org/。
另外,我还成功地使用了:ACADO,http://sourceforge.net/p/acado/wiki/Home/(在C++中)
3ks5zfa02#
现在有一些非常强大的建模工具,如CasADi、Pyomo和Gekko,这些工具在提出这个问题时还不存在。这里有一个Gekko的解决方案。一个问题是
sqrt(u)
需要有一个正的u
值,以避免虚数。