将 GEKKO ```Var``` 的初始值和最终值限制为基于数据的曲线
constrain the initial and final values of a GEKKO ```Var``` to a data-based curve
我正在尝试解决地球轨道的低推力优化控制问题，即从一个轨道到另一个轨道。问题的表述包括六个状态（r_1、r_2、r_3、v_1、v_2、v_3）和3个控制（u_1, u_2, u_3) 使用简化的重力点模型。当我指定完整的初始状态和一半的最终状态时，求解器收敛并产生一个好的解决方案。当我尝试完整的最终状态时，问题是过度约束。
我的补救方法是让轨迹在轨道曲线上的任意点离开初始轨道，并在最终轨道曲线上的任意点加入最终轨道，赋予它更多的自由度。 有没有办法将所有 6 个状态的初始值和最终值限制为 cspline 曲线？ 这是我目前尝试过的方法：
m = GEKKO(remote=True)
m.time = np.linspace(0, tof, t_steps)
theta_i = m.FV(lb = 0, ub = 360)
theta_f = m.FV(lb = 0, ub = 360)
rx_i = m.Param()
ry_i = m.Param()
rz_i = m.Param()
vx_i = m.Param()
vy_i = m.Param()
vz_i = m.Param()
rx_f = m.Param()
ry_f = m.Param()
rz_f = m.Param()
vx_f = m.Param()
vy_f = m.Param()
vz_f = m.Param()
m.cspline(theta_i, rx_i, initial.theta, initial.r[:, 0])
m.cspline(theta_i, ry_i, initial.theta, initial.r[:, 1])
m.cspline(theta_i, rz_i, initial.theta, initial.r[:, 2])
m.cspline(theta_i, vx_i, initial.theta, initial.v[:, 0])
m.cspline(theta_i, vy_i, initial.theta, initial.v[:, 1])
m.cspline(theta_i, vz_i, initial.theta, initial.v[:, 2])
m.cspline(theta_f, rx_f, final.theta, final.r[:, 0])
m.cspline(theta_f, ry_f, final.theta, final.r[:, 1])
m.cspline(theta_f, rz_f, final.theta, final.r[:, 2])
m.cspline(theta_f, vx_f, final.theta, final.v[:, 0])
m.cspline(theta_f, vy_f, final.theta, final.v[:, 1])
m.cspline(theta_f, vz_f, final.theta, final.v[:, 2])
r1 = m.Var(rx_i)
r2 = m.Var(ry_i)
r3 = m.Var(rz_i)
r1dot = m.Var(vx_i)
r2dot = m.Var(vy_i)
r3dot = m.Var(vz_i)
u1 = m.Var(lb = -max_u, ub = max_u)
u2 = m.Var(lb = -max_u, ub = max_u)
u3 = m.Var(lb = -max_u, ub = max_u)
m.Equation(r1.dt() == r1dot)
m.Equation(r2.dt() == r2dot)
m.Equation(r3.dt() == r3dot)
r = m.Intermediate(m.sqrt(r1**2 + r2**2 + r3**2))
v = m.Intermediate(m.sqrt(r1dot**2 + r2dot**2 + r3dot**3))
m.Equation(-mu*r1/r**3 == r1dot.dt() + u1)
m.Equation(-mu*r2/r**3 == r2dot.dt() + u2)
m.Equation(-mu*r3/r**3 == r3dot.dt() + u3)
m.fix_final(r1, rx_f)
m.fix_final(r2, ry_f)
m.fix_final(r3, rz_f)
m.fix_final(r1dot, vx_f)
m.fix_final(r2dot, vy_f)
m.fix_final(r3dot, vz_f)
m.Minimize(m.integral(u1**2 + u2**2 + u3**2))
m.options.IMODE = 6
m.options.solver = 3
#m.options.ATOL = 1e-3
m.options.MAX_ITER = 300
m.solve(disp=True) # solve
使用 cspline，我试图让 GEKKO 从我生成的关于采样真实异常状态的数据中选择真实异常的固定值（该参数决定了你在轨道上的距离）会有相关的位置和速度状态。任何帮助将不胜感激！
解决方案：
我将最终约束实现为“软约束”，即要最小化的数量，而不是使用 fix_final 的硬约束。具体实现如下
final = np.zeros(len(m.time))
final[-1] = 1
final = m.Param(value=final)
m.Obj(final*(r1-final_o.r[0, 0])**2)
m.Obj(final*(r2-final_o.r[0, 1])**2)
m.Obj(final*(r3-final_o.r[0, 2])**2)
m.Obj(final*(r1dot-final_o.v[0, 0])**2)
m.Obj(final*(r2dot-final_o.v[0, 1])**2)
m.Obj(final*(r3dot-final_o.v[0, 2])**2)
final 使求解器在相乘时只考虑最后一项（结束边界）。
优化器通常更难准确地达到固定端点，尤其是当它依赖于复杂的移动序列时。这通常会导致不可行的解决方案。另一种方法是创建一个软约束（objective 最小化）来惩罚与最终轨迹的偏差。这是一个类似的例子：
import matplotlib.animation as animation
import numpy as np
from gekko import GEKKO
#Defining a model
m = GEKKO()
#################################
#Weight of item
m2 = 1
#################################
#Defining the time, we will go beyond the 6.2s
#to check if the objective was achieved
m.time = np.linspace(0,8,100)
end_loc = int(100.0*6.2/8.0)
#Parameters
m1a = m.Param(value=10)
m2a = m.Param(value=m2)
final = np.zeros(len(m.time))
for i in range(len(m.time)):
if m.time[i] < 6.2:
final[i] = 0
else:
final[i] = 1
final = m.Param(value=final)
#MV
ua = m.Var(value=0)
#State Variables
theta_a = m.Var(value=0)
qa = m.Var(value=0)
ya = m.Var(value=-1)
va = m.Var(value=0)
#Intermediates
epsilon = m.Intermediate(m2a/(m1a+m2a))
#Defining the State Space Model
m.Equation(ya.dt() == va)
m.Equation(va.dt() == -epsilon*theta_a + ua)
m.Equation(theta_a.dt() == qa)
m.Equation(qa.dt() == theta_a -ua)
#Define the Objectives
#Make all the state variables be zero at time >= 6.2
m.Obj(final*ya**2)
m.Obj(final*va**2)
m.Obj(final*theta_a**2)
m.Obj(final*qa**2)
m.fix(ya,pos=end_loc,val=0.0)
m.fix(va,pos=end_loc,val=0.0)
m.fix(theta_a,pos=end_loc,val=0.0)
m.fix(qa,pos=end_loc,val=0.0)
#Try to minimize change of MV over all horizon
m.Obj(0.001*ua**2)
m.options.IMODE = 6 #MPC
m.solve() #(disp=False)
这个 example 使用软约束和硬约束的组合来帮助它找到最佳解决方案。
我正在尝试解决地球轨道的低推力优化控制问题，即从一个轨道到另一个轨道。问题的表述包括六个状态（r_1、r_2、r_3、v_1、v_2、v_3）和3个控制（u_1, u_2, u_3) 使用简化的重力点模型。当我指定完整的初始状态和一半的最终状态时，求解器收敛并产生一个好的解决方案。当我尝试完整的最终状态时，问题是过度约束。
我的补救方法是让轨迹在轨道曲线上的任意点离开初始轨道，并在最终轨道曲线上的任意点加入最终轨道，赋予它更多的自由度。 有没有办法将所有 6 个状态的初始值和最终值限制为 cspline 曲线？ 这是我目前尝试过的方法：
m = GEKKO(remote=True)
m.time = np.linspace(0, tof, t_steps)
theta_i = m.FV(lb = 0, ub = 360)
theta_f = m.FV(lb = 0, ub = 360)
rx_i = m.Param()
ry_i = m.Param()
rz_i = m.Param()
vx_i = m.Param()
vy_i = m.Param()
vz_i = m.Param()
rx_f = m.Param()
ry_f = m.Param()
rz_f = m.Param()
vx_f = m.Param()
vy_f = m.Param()
vz_f = m.Param()
m.cspline(theta_i, rx_i, initial.theta, initial.r[:, 0])
m.cspline(theta_i, ry_i, initial.theta, initial.r[:, 1])
m.cspline(theta_i, rz_i, initial.theta, initial.r[:, 2])
m.cspline(theta_i, vx_i, initial.theta, initial.v[:, 0])
m.cspline(theta_i, vy_i, initial.theta, initial.v[:, 1])
m.cspline(theta_i, vz_i, initial.theta, initial.v[:, 2])
m.cspline(theta_f, rx_f, final.theta, final.r[:, 0])
m.cspline(theta_f, ry_f, final.theta, final.r[:, 1])
m.cspline(theta_f, rz_f, final.theta, final.r[:, 2])
m.cspline(theta_f, vx_f, final.theta, final.v[:, 0])
m.cspline(theta_f, vy_f, final.theta, final.v[:, 1])
m.cspline(theta_f, vz_f, final.theta, final.v[:, 2])
r1 = m.Var(rx_i)
r2 = m.Var(ry_i)
r3 = m.Var(rz_i)
r1dot = m.Var(vx_i)
r2dot = m.Var(vy_i)
r3dot = m.Var(vz_i)
u1 = m.Var(lb = -max_u, ub = max_u)
u2 = m.Var(lb = -max_u, ub = max_u)
u3 = m.Var(lb = -max_u, ub = max_u)
m.Equation(r1.dt() == r1dot)
m.Equation(r2.dt() == r2dot)
m.Equation(r3.dt() == r3dot)
r = m.Intermediate(m.sqrt(r1**2 + r2**2 + r3**2))
v = m.Intermediate(m.sqrt(r1dot**2 + r2dot**2 + r3dot**3))
m.Equation(-mu*r1/r**3 == r1dot.dt() + u1)
m.Equation(-mu*r2/r**3 == r2dot.dt() + u2)
m.Equation(-mu*r3/r**3 == r3dot.dt() + u3)
m.fix_final(r1, rx_f)
m.fix_final(r2, ry_f)
m.fix_final(r3, rz_f)
m.fix_final(r1dot, vx_f)
m.fix_final(r2dot, vy_f)
m.fix_final(r3dot, vz_f)
m.Minimize(m.integral(u1**2 + u2**2 + u3**2))
m.options.IMODE = 6
m.options.solver = 3
#m.options.ATOL = 1e-3
m.options.MAX_ITER = 300
m.solve(disp=True) # solve
使用 cspline，我试图让 GEKKO 从我生成的关于采样真实异常状态的数据中选择真实异常的固定值（该参数决定了你在轨道上的距离）会有相关的位置和速度状态。任何帮助将不胜感激！
解决方案：
我将最终约束实现为“软约束”，即要最小化的数量，而不是使用 fix_final 的硬约束。具体实现如下
final = np.zeros(len(m.time))
final[-1] = 1
final = m.Param(value=final)
m.Obj(final*(r1-final_o.r[0, 0])**2)
m.Obj(final*(r2-final_o.r[0, 1])**2)
m.Obj(final*(r3-final_o.r[0, 2])**2)
m.Obj(final*(r1dot-final_o.v[0, 0])**2)
m.Obj(final*(r2dot-final_o.v[0, 1])**2)
m.Obj(final*(r3dot-final_o.v[0, 2])**2)
final 使求解器在相乘时只考虑最后一项（结束边界）。
优化器通常更难准确地达到固定端点，尤其是当它依赖于复杂的移动序列时。这通常会导致不可行的解决方案。另一种方法是创建一个软约束（objective 最小化）来惩罚与最终轨迹的偏差。这是一个类似的例子：
import matplotlib.animation as animation
import numpy as np
from gekko import GEKKO
#Defining a model
m = GEKKO()
#################################
#Weight of item
m2 = 1
#################################
#Defining the time, we will go beyond the 6.2s
#to check if the objective was achieved
m.time = np.linspace(0,8,100)
end_loc = int(100.0*6.2/8.0)
#Parameters
m1a = m.Param(value=10)
m2a = m.Param(value=m2)
final = np.zeros(len(m.time))
for i in range(len(m.time)):
if m.time[i] < 6.2:
final[i] = 0
else:
final[i] = 1
final = m.Param(value=final)
#MV
ua = m.Var(value=0)
#State Variables
theta_a = m.Var(value=0)
qa = m.Var(value=0)
ya = m.Var(value=-1)
va = m.Var(value=0)
#Intermediates
epsilon = m.Intermediate(m2a/(m1a+m2a))
#Defining the State Space Model
m.Equation(ya.dt() == va)
m.Equation(va.dt() == -epsilon*theta_a + ua)
m.Equation(theta_a.dt() == qa)
m.Equation(qa.dt() == theta_a -ua)
#Define the Objectives
#Make all the state variables be zero at time >= 6.2
m.Obj(final*ya**2)
m.Obj(final*va**2)
m.Obj(final*theta_a**2)
m.Obj(final*qa**2)
m.fix(ya,pos=end_loc,val=0.0)
m.fix(va,pos=end_loc,val=0.0)
m.fix(theta_a,pos=end_loc,val=0.0)
m.fix(qa,pos=end_loc,val=0.0)
#Try to minimize change of MV over all horizon
m.Obj(0.001*ua**2)
m.options.IMODE = 6 #MPC
m.solve() #(disp=False)
这个 example 使用软约束和硬约束的组合来帮助它找到最佳解决方案。