如何以图形方式验证我的动态优化结果，与 gekko 中的初始条件进行比较
Posted
技术标签:
【中文标题】如何以图形方式验证我的动态优化结果，与 gekko 中的初始条件进行比较【英文标题】：How to verify my dynamic optimization results graphically, comparing with initial conditions in gekko 【发布时间】：2021-11-17 01:25:50 【问题描述】：朋友们，Hedengren 教授早上好，我是 Python 新手，对 Gekko 更是新手，首先我想知道我在 Gekko 中的代码是否正确，即根据我的身体需要，考虑到我的方程式是正确的。
我的模型尝试优化变量M2 和l_M2（或这两个变量的组合），以便在模块中最小化我的变量q1 的运动幅度（正或负），我的模型从放置的文本文件here 接收输入，模型解决方案必须遵守以下内容：
M2和l_M2的初始值，求解模型，得到q1的最大幅度（正或负）；
输入值不随水平变化；
在每次迭代中，变量c_m2 的值必须根据M2 和l_M2 的值进行更新，并且必须在整个范围内保持不变。
为了最小化变量q1我提出了两种类型的目标，我不会同时使用它们：
1000*q1**2;
最小化x1 = integral (0.5 q1 ** 2)dt从0到t的积分，为其创建一个辅助变量x1。
有待解决的疑问
求解模型时发现c_m2的值（在初始点）为0。这是不正确的，因为它应该与下面的值相同，所以我的代码有错误，我不知道。如何解决；
另一方面，我希望能够将模型的响应与变量的初始值与优化值的响应进行比较（如图所示），但我不能了解如何使用初始值保存我的响应。
Optimization check figure
在这种情况下使用m.options.IMODE = 6 是否正确？
这是我的代码：
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
###################### CREATION OF LOAD RECORD
filename= 'Inputs 0.02sec.txt'
input_l=(np.loadtxt(filename, skiprows=1, dtype=float).flatten()).tolist()
dt=0.02
len_inputs=len(input_l)
m=GEKKO()
# time vector
t_final=dt*(len_inputs-1)
m.time=np.linspace(0, t_final, len_inputs)
# parameters
M1=m.Param(value=21956548.3771968)
Ri=m.Param(value=10609404.1758615)
taxa1=m.Param(value=0.02)
taxa2=m.Param(value=0.005)
grv=m.Param(value=9.80665)
in_loads=m.Param(value=input_l)
m.options.NODES = 4
m.options.IMODE = 6 #MPC
#Intermedias
Om1=m.Intermediate(m.sqrt(Ri/M1))
C_M1=m.Intermediate(2*M1*Om1*taxa1)
# variables
M2=m.FV(value=0.10*21956548.3771968,lb=0.01*M1 , ub=0.20*M1)
M2.STATUS = 1
l_M2=m.FV(value=7, lb=1, ub=20)
l_M2.STATUS = 1
c_m2=m.Var(value=2*taxa2*M2*m.sqrt(grv/l_M2))
x1=m.Var(value=0) # auxiliar variable for integral of x1=0.5*integral(q1**2)dt
q1=m.Var(value=0)
q1_p=m.Var(value=0)
q2=m.Var(value=0)
q2_p=m.Var(value=0)
# auxiliar equation for minimization of integral of x1=0.5*integral(q1**2)dt
m.Equation(x1.dt()==0.5*(q1**2))
# equations for actualization of c_m2
m.Equation(c_m2==2*taxa2*m.sqrt(grv/l_M2))
# equations of state
m.Equation(q1.dt()==q1_p)
m.Equation(q1_p.dt()==((-Ri*q1-C_M1*q1_p+M2*grv*q2+(c_m2*q2_p)/l_M2) \
/M1-in_loads))
m.Equation(q2.dt()==q2_p)
m.Equation(q2_p.dt()==(Ri*q1+C_M1*q1_p-(M1+M2)*grv*q2)/(l_M2*M1) \
-c_m2*(M1+M2)*q2_p/(M1*M2*l_M2**2))
m.Obj(1000*q1**2) # for minimization of q1 (1000*q1**2)
# m.Obj(x1) # for minimization of integral 0.5*q1**2
m.solve()
######################################### Plotting the results
fig=plt.figure(1)
ax4 = fig.add_subplot(1,1,1)
ax4.plot(m.time, q1.value, ls='-', label=f'q1 Opt')
ax4.set_ylabel('Amplitude of q1 [m]')
ax4.set_xlabel('Time [sec]')
ax4.set_title('Time - Amplitude \n')
ax4.legend(loc='best')
plt.grid()
minimo,maximo=min(q1.value),max(q1.value)
Max_q1=max(abs(minimo),abs(maximo))
# print results
print ('')
print ('--- Results of the Optimization Problem ---')
print ('M2= ' + str(M2.value))
print ('l_M2 = ' + str(l_M2.value))
print ('c_m2 = ' + str(c_m2.value))
print ('Absolute Max Amplitude q1= ', Max_q1)
print ('Percentage of massa m2= ' + str(M2.value[-1]/M1.value[-1]))
plt.show()
【问题讨论】：
有m.integral(0.5 q1**2) 功能，如果使用它来定义您的问题更容易。另外，请尝试使用m.Minimize() 而不是m.Obj() 以提高可读性。
【参考方案1】：
在这个应用程序上做得很好。你描述的一切看起来都是正确的。 c_m2 的值在开始时为零，因为初始条件固定在指定值（或默认值 0），并且在 t=0 处无法求解方程。
尝试使用两个m.solve() 命令，第一个使用m.options.COLDSTART=1（临时设置STATUS=0）来查看初始和优化的解决方案。
m.options.COLDSTART = 1
m.solve()
fig=plt.figure(1)
ax4 = fig.add_subplot(1,1,1)
ax4.plot(m.time, q1.value, ls='-', label=f'q1 Initial')
m.options.COLDSTART = 0
m.options.TIME_SHIFT = 0
m.solve()
ax4.plot(m.time, q1.value, ls='-', label=f'q1 Opt')
ax4.set_ylabel('Amplitude of q1 [m]')
ax4.set_xlabel('Time [sec]')
ax4.set_title('Time - Amplitude \n')
ax4.legend(loc='best')
这是完整的脚本：
import numpy as np
import matplotlib.pyplot as plt
from gekko import GEKKO
###################### CREATION OF LOAD RECORD
filename= 'Inputs 0.02sec.txt'
input_l=(np.loadtxt(filename, skiprows=1, dtype=float).flatten()).tolist()
dt=0.02
len_inputs=len(input_l)
m=GEKKO()
# time vector
t_final=dt*(len_inputs-1)
m.time=np.linspace(0, t_final, len_inputs)
# parameters
M1=m.Param(value=21956548.3771968)
Ri=m.Param(value=10609404.1758615)
taxa1=m.Param(value=0.02)
taxa2=m.Param(value=0.005)
grv=m.Param(value=9.80665)
in_loads=m.Param(value=input_l)
m.options.NODES = 4
m.options.IMODE = 6 #MPC
#Intermedias
Om1=m.Intermediate(m.sqrt(Ri/M1))
C_M1=m.Intermediate(2*M1*Om1*taxa1)
# variables
M2=m.FV(value=0.10*21956548.3771968,lb=0.01*M1 , ub=0.20*M1)
M2.STATUS = 1
l_M2=m.FV(value=7, lb=1, ub=20)
l_M2.STATUS = 1
c_m2=m.Var(value=2*taxa2*M2*m.sqrt(grv/l_M2))
x1=m.Var(value=0) # auxiliar variable for integral of x1=0.5*integral(q1**2)dt
q1=m.Var(value=0)
q1_p=m.Var(value=0)
q2=m.Var(value=0)
q2_p=m.Var(value=0)
# auxiliar equation for minimization of integral of x1=0.5*integral(q1**2)dt
m.Equation(x1.dt()==0.5*(q1**2))
# equations for actualization of c_m2
m.Equation(c_m2==2*taxa2*m.sqrt(grv/l_M2))
# equations of state
m.Equation(q1.dt()==q1_p)
m.Equation(q1_p.dt()==((-Ri*q1-C_M1*q1_p+M2*grv*q2+(c_m2*q2_p)/l_M2) \
/M1-in_loads))
m.Equation(q2.dt()==q2_p)
m.Equation(q2_p.dt()==(Ri*q1+C_M1*q1_p-(M1+M2)*grv*q2)/(l_M2*M1) \
-c_m2*(M1+M2)*q2_p/(M1*M2*l_M2**2))
m.Minimize(1000*q1**2) # for minimization of q1 (1000*q1**2)
m.options.COLDSTART = 1
m.solve()
fig=plt.figure(1)
ax4 = fig.add_subplot(1,1,1)
ax4.plot(m.time, q1.value, ls='-', label=f'q1 Initial')
m.options.COLDSTART = 0
m.options.TIME_SHIFT = 0
m.solve()
ax4.plot(m.time, q1.value, ls='-', label=f'q1 Opt')
ax4.set_ylabel('Amplitude of q1 [m]')
ax4.set_xlabel('Time [sec]')
ax4.set_title('Time - Amplitude \n')
ax4.legend(loc='best')
minimo,maximo=min(q1.value),max(q1.value)
Max_q1=max(abs(minimo),abs(maximo))
# print results
print ('')
print ('--- Results of the Optimization Problem ---')
print ('M2= ' + str(M2.value))
print ('l_M2 = ' + str(l_M2.value))
print ('c_m2 = ' + str(c_m2.value))
print ('Absolute Max Amplitude q1= ', Max_q1)
print ('Percentage of massa m2= ' + str(M2.value[-1]/M1.value[-1]))
plt.show()
【讨论】：
教授，请教一个问题，我有一个已知的最佳点，当我将其设置为变量M2 = m.FV (value = 0.20 * 21956548.3771968, lb = 0.01 * M1, ub = 0.20 * M1)和l_M2 = m .FV (value = 13.2927529387356, lb = 1, ub = 20)的起点时，结果比已知的最佳点差，为什么会这样？有什么办法可以防止我得到这个错误的结果？还是将同样的最佳点还给我？最大结果是：绝对最大振幅 q1 Opt = 0.22987254086 绝对最大振幅 q1 初始条件 = 0.20725808504
另一方面，当我从接近这个已知点的初始条件开始时，我无法达到这个最佳点，因此我的结果更糟。
你想最大化 q1 吗？尝试使用m.maximize() 而不是最小化的m.Obj()。还要检查初始猜测时是否满足约束。
没有教授，我正在尝试最小化q1的最大幅度，初始条件满足约束，但是求解时，优化了q1的最大值大于初始条件的最大值，正如我在第一条评论中提到的，我唯一做的就是修改了两个变量的初始条件您输入的代码 M2 = m.FV (value = 0.20 * 21956548.3771968, lb = 0.01 * M1, ub = 0.20 * M1) 和 l_M2 = m .FV (value = 13.2927529387356, lb = 1, ub = 20) 。提前教授，非常感谢您的帮助。
M2的初始值与其上限相同。以上是关于如何以图形方式验证我的动态优化结果，与 gekko 中的初始条件进行比较的主要内容，如果未能解决你的问题，请参考以下文章
GEKKO - 如何修复 Python Gekko Max Equation 错误 - 元素数