Jacob Olsen
Reputation: 103
Using Gekko to reproduce a control history in an observed trajectory
I am trying to get Gekko to reproduce the bank angle history from the Apollo 10 reentry event. I have altitude, velocity, and latitude state histories and the actual NASA bank angle data for comparison. My path constraints look like this:
m.Minimize(1e3*(r_desired-r)**2)
m.Minimize(1e6*(V_desired-V)**2)
m.Minimize((lat_desired-phi)**2)
m.Minimize(0.5*bank**2)
Couple of questions:
Did I implement the path constraint correctly? Is there a smarter way to do that?
I am getting a converged solution but for the life of me I cannot get the bank angle to match. It gets pretty close at times (especially if you consider 180=-180) but definitely off at other times. Any recommendations on how to get better agreement?
I have tried tons of different combinations and iterations to include weighting the constraints, using linear interpolation vs using a polynomial fit for the historical data, turning on and off other constraints, limiting the bounds on the bank angle, etc. Any help would be very appreciated. Thanks!
Code and state history data are below.
from gekko import GEKKO
import numpy as np
import matplotlib.pyplot as plt
import math
import pandas as pd
pi = math.pi
#======================================#
#OPTIONS
#======================================#
rotating_earth = 0 #### 0 = 'no' 1 = 'yes'####
#======================================#
#Gekko setup
#======================================#
m = GEKKO() # initialize GEKKO
tfin = 200
nt = 201 #simulation time is 2500 seconds
gekkoTime = np.linspace(0,tfin,nt) #time array
m.time = gekkoTime
#======================================#
#Params
#======================================#
r_earth = m.Param(value=6378.137) #radius of the earth
g0 = m.Param(value=9.80665/1000) #acceleration of gravity
ro_s = m.Param(value=1.225e9) #atmospheric density at sea level
beta = m.Param(value=0.14) #inverse scale height
omega = m.Param(value=0) #rotation of the earth, rad/sec
mass = m.Param(value=5498.2) # mass at reentry
S = m.Param(value=12.02e-6) #reference area km^2
Cl = m.Param(value=0.4082) #coefficient of lift
Cd = m.Param(value=1.2569) #drag
CS = m.Param(value=Cl) #coeff of sideslip
T = m.Param(0) #N, thrust of satellite
pcsi = m.Param(0) #rad, "Vertical" thrust vector angle
epsilon = m.Param(0) #rad, "Horizontal" thrust vector angle
#======================================#
#Initial Conditions
#======================================#
#States at entry
Re = 6498.270 #km
flp_e= -6.6198381*(pi/180) #flight path angle at entry
psi_e = 18.0683*(pi/180) #heading angle
V_e = 11.06715 #velocity km/sec
theta_e = 174.24384*(pi/180) #longitude east
phi_e = -23.51457*(pi/180) #latitude south
if rotating_earth == 1: # if the rotating Earth is turned on
V_e = m.Param(10.6589) #entry Velocity, km/sec
flp_e = m.Param(-6.87459*pi/180) #initial flight path angle
psi_e = m.Param(19.9425*pi/180) #initial heading angle
Azi_e = m.Param(psi_e+(pi/2-2*psi_e)) #intial azimuth angle
omega = m.Param(value=7.2921159e-5) #rotation of the earth, rad/sec
#======================================#
#Variables
#======================================#
# #state variables
# r = m.Var(value=Re, lb=6378.137, ub=Re)
# V = m.Var(value=V_e, lb=0, ub=V_e)
# phi = m.Var(value=phi_e, lb=-90*pi/180, ub=70*pi/180) #latitude
# theta = m.Var(value=theta_e, lb=-pi, ub=pi) #longitude
# flp = m.Var(value=flp_e, lb=-90*pi/180, ub=10*pi/180)
# psi = m.Var(value=psi_e)
#state variables
# r = m.Var(value=Re, lb=6378.137, ub=Re)
# V = m.Var(value=V_e, lb=0, ub=V_e)
# phi = m.Var(value=phi_e, lb=-25*pi/180, ub=-13*pi/180) #latitude
# theta = m.Var(value=theta_e, lb=150*pi/180, ub=220*pi/180) #longitude
# flp = m.Var(value=flp_e, lb=-90*pi/180, ub=10*pi/180)
# psi = m.Var(value=psi_e)
#state variables
r = m.Var(value=Re)
V = m.Var(value=V_e)
phi = m.Var(value=phi_e) #latitude
theta = m.Var(value=theta_e) #longitude
flp = m.Var(value=flp_e)
psi = m.Var(value=psi_e)
#control variables
bank = m.MV(value=0, lb=-180*pi/180, ub=pi) # bank angle, rad
# bank = m.MV() # bank angle, rad
bank.STATUS = 1
#===============================#
#APOLLO HISTORICAL DATA
#===============================#
apolloAlt = np.array(pd.read_excel('ApolloData.xls', sheet_name="Altitude"))
apolloVel = np.array(pd.read_excel('ApolloData.xls', sheet_name="Velocity"))
apolloLat = np.array(pd.read_excel('ApolloData.xls', sheet_name="Latitude"))
apolloBank = np.array(pd.read_excel('ApolloData.xls', sheet_name="Bank"))
def lin_interp(gekkoTime, ApolloData):
return np.array(np.interp(gekkoTime, ApolloData[:,0], ApolloData[:,1]))
# alt_desired = lin_interp(gekkoTime, apolloAlt)
# r_desired = alt_desired+r_earth.value
Vel_desired = lin_interp(gekkoTime,apolloVel)
# lat_desired = lin_interp(gekkoTime,apolloLat)
bankTime = np.linspace(0,500,2000)
bank_desired = lin_interp(bankTime,apolloBank)
#curve fit
x1 = apolloAlt[:,0]
y2= apolloAlt[:,1]
p30 = np.poly1d(np.polyfit(x1, y2, 30))
alt_desired = p30(gekkoTime)
r_desired = alt_desired+r_earth.value
# x1 = apolloVel[:,0]
# y2 = apolloVel[:,1]
# p30 = np.poly1d(np.polyfit(x1, y2, 30))
# # Vel_desired = p30(gekkoTime)
x1 = apolloLat[:,0]
y2 = apolloLat[:,1]
p30 = np.poly1d(np.polyfit(x1, y2, 30))
lat_desired = p30(gekkoTime)
path = m.Param(value=alt_desired)
#======================================#
#Intermediate Variables
#======================================#
g = m.Intermediate(g0*(r_earth/r)**2) #calculate gravity
ro = m.Intermediate(ro_s*m.exp(-beta*(r-r_earth))) #calculate ro
D = m.Intermediate(((ro*Cd*S)/2)*V**2) #calculate drag
L = m.Intermediate(((ro*Cl*S)/2)*V**2) #calculate lift
# error = m.Intermediate(m.abs3(path-r))
#======================================#
#EQUATIONS OF MOTION
#======================================#
m.Equation(r.dt() == V*m.sin(flp)) #calculate r.dt()
m.Equation(theta.dt() == (V*m.cos(flp)*m.cos(psi))/(r*m.cos(phi))) #calculate theta .dt()
m.Equation(phi.dt() == (V*m.cos(flp)*m.sin(psi))/r) #calculate phi .dt()
m.Equation(V.dt() == (T/mass)*(m.cos(pcsi)*m.cos(epsilon))-(D/mass) - g*m.sin(flp)+r*omega**2*m.cos(phi)\
*(m.cos(phi)*m.sin(flp)-m.sin(phi)*m.sin(psi)*m.cos(flp))) #calculate V.dt()
m.Equation(flp.dt() == ((T/mass)*(m.sin(pcsi)*m.sin(bank)+m.cos(pcsi)*m.sin(epsilon)*m.cos(bank))+(L/mass)\
*m.cos(bank)-g*m.cos(flp)+((V**2)/r)*m.cos(flp)+2*V*omega*m.cos(phi)*m.cos(psi)+r*omega**2*m.cos(phi)\
*(m.cos(phi)*m.cos(flp)+m.sin(phi)*m.sin(psi)*m.sin(flp)))/V) #calculate flp .dt()
m.Equation(psi.dt() == ((1/(mass*m.cos(flp)))*(T*(m.cos(pcsi)*m.sin(epsilon)*m.sin(bank)-m.sin(pcsi)\
*m.cos(bank))+L*m.sin(bank))-((V**2)/r)*m.cos(flp)*m.cos(psi)*m.tan(phi)+2*V*omega\
*(m.sin(psi)*m.cos(phi)*m.tan(flp)-m.sin(phi))-((r*omega**2)/m.cos(flp))*m.sin(phi)\
*m.cos(phi)*m.cos(psi))/V) #calculate psi .dt()
# m.Equation(r.dt() == V*m.sin(flp)) #calculate r.dt()
#
# m.Equation(theta.dt()*(r*m.cos(phi)) == (V*m.cos(flp)*m.cos(psi))) #calculate theta .dt()
#
# m.Equation(phi.dt()*(r) == (V*m.cos(flp)*m.sin(psi))) #calculate phi .dt()
#
# m.Equation(V.dt() == -D/mass-g*m.sin(flp)) #calculate V.dt()
#
# m.Equation(flp.dt()*V == ((L/mass)*m.cos(bank) - g*m.cos(flp)+((V**2)/r)*m.cos(flp))) #calculate gamma .dt()
#
# m.Equation(psi.dt()*V == ((L*m.sin(bank))/(mass*m.cos(flp))-((V**2)/r)*m.cos(flp)*m.cos(psi)*m.tan(phi))) #calculate psi .dt()
#======================================#
#Optimize
#======================================#
p = np.zeros(nt) # mark final time point
p[-1] = 1.0
final = m.Param(value=p)
r_desired= m.Param(value=r_desired)
V_desired= m.Param(value=Vel_desired)
lat_desired = m.Param(value=lat_desired*pi/180)
# m.Minimize(final*r)
# m.Minimize(error)
# m.Minimize(.5*bank**2)
# m.Minimize(m.abs(r_desired-r))
m.Minimize(1e3*(r_desired-r)**2)
m.Minimize(1e6*(V_desired-V)**2)
m.Minimize((lat_desired-phi)**2)
m.Minimize(0.5*bank**2)
# m.Minimize(m.abs(bank))
# m.options.MAX_ITER=1000
m.options.IMODE = 6
m.options.SOLVER = 3
m.solve()
#======================================#
#Plotting
#======================================#
alt = np.array(r.value)
alt = alt-r_earth.value
vel = np.array(V.value)
latitude = np.array(phi.value)*180/pi
longitude = np.array(theta.value)*180/pi
fpa = np.array(flp.value)*180/pi
psi = np.array(psi.value)*180/pi
bank = np.array(bank.value)*180/pi
fig = plt.figure(figsize=(12,8))
ax1 = fig.add_subplot(421)
ax1.plot(m.time,alt,label='Alt')
ax1.plot(m.time,alt_desired,label='Nasa Alt')
# ax1.title.set_text('Altitude (km)')
ax1.set_xlabel('Time (sec)')
ax1.set_ylabel('Altitude (km)')
ax1.grid(True)
plt.legend()
ax2 = fig.add_subplot(422)
ax2.plot(m.time,vel,label='Vel')
ax2.plot(m.time,V_desired,label='Nasa Vel')
ax2.set_xlabel('Time (sec)')
ax2.set_ylabel('Vel (km/sec)')
ax2.grid(True)
plt.legend()
ax3 = fig.add_subplot(423)
ax3.plot(m.time,latitude)
ax3.set_xlabel('Time (sec)')
ax3.set_ylabel('Latitude (deg)')
ax3.grid(True)
ax4 = fig.add_subplot(424)
ax4.plot(m.time,longitude)
ax4.set_xlabel('Time (sec)')
ax4.set_ylabel('Longitude (deg)')
ax4.grid(True)
ax5 = fig.add_subplot(425)
ax5.plot(m.time,fpa)
ax5.set_xlabel('Time (sec)')
ax5.set_ylabel('FPA (deg)')
ax5.grid(True)
ax6 = fig.add_subplot(426)
ax6.plot(m.time,psi)
ax6.set_xlabel('Time (sec)')
ax6.set_ylabel('Azimuth')
ax6.grid(True)
ax8 = fig.add_subplot(427)
ax8.plot(m.time, bank, label='Bank')
ax8.plot(bankTime, bank_desired,label='Nasa Bank')
ax8.set_xlabel('Time (sec)')
ax8.set_ylabel('Bank (Deg)')
ax8.grid(True)
plt.legend(loc='upper right')
plt.suptitle("Apollo 10 Trajectory Matching")
plt.savefig('Apollo10_TrajMatching.png')
plt.show()
# plt.subplot(4,2,1)
# plt.plot(m.time,alt,label='alt')
# plt.plot(m.time,alt_desired, label='Nasa Alt')
# plt.legend()
# plt.subplot(4,2,3)
# plt.plot(m.time,latitude,label='lat')
# plt.legend()
# plt.subplot(4,2,4)
# plt.plot(m.time,longitude,label='lon')
# plt.legend()
# plt.subplot(4,2,2)
# plt.plot(m.time,vel,label='vel')
# plt.plot(m.time,V_desired,label='Nasa Vel')
# plt.legend()
# plt.subplot(4,2,5)
# plt.plot(m.time,fpa,label='fpa')
# plt.legend()
# plt.subplot(4,2,6)
# plt.plot(m.time,psi,label='heading')
# plt.legend()
# plt.subplot(4,2,7)
# plt.plot(m.time,bank,label='bank')
# plt.plot(m.time,bank_desired,label='Nasa Bank')
# plt.xlabel('Time')
# plt.legend()
# plt.show()
and State history:
Time Altitude Time Velocity Time Latitude Time Bank Angle
0 123.936 11.8033 11.0962 0.101882 -23.4628 0 0
1.85811 118.083 19.6721 11.0962 9.63964 -23.1728 1 0
5.57432 112.761 33.4426 11.0962 19.1774 -22.8829 2 0
9.29054 108.502 41.3115 11.046 30.6196 -22.5514 3 0
13.0068 103.711 47.2131 10.9958 43.9818 -22.0959 4 0
16.723 99.4519 53.1148 10.8954 55.4318 -21.7231 5 0
18.5811 95.7267 57.0492 10.7448 64.9853 -21.3505 6 0
22.2973 92.5316 62.9508 10.5941 78.3318 -20.9776 7 0
26.0135 88.8046 66.8852 10.4435 85.973 -20.6878 8 0
29.7297 85.6095 70.8197 10.1925 95.4951 -20.4804 9 0
31.5878 82.4163 72.7869 9.94142 105.025 -20.2318 10 0
35.3041 80.285 74.7541 9.74059 112.651 -20.0246 11 0
37.1622 76.5598 78.6885 9.48954 120.284 -19.7761 12 0
40.8784 73.3647 84.5902 9.23849 127.902 -19.6102 13 0
46.4527 70.6998 86.5574 9.08787 135.527 -19.403 14 0
50.1689 67.5047 92.459 8.78661 141.248 -19.2373 15 0
53.8851 64.8415 98.3607 8.53556 148.866 -19.0715 16 0
57.6014 62.1783 104.262 8.28452 154.587 -18.9057 17 0
61.3176 60.0471 112.131 8.03347 160.316 -18.6987 18 0
63.1757 58.9815 118.033 7.88285 173.639 -18.4498 19 0
65.0338 57.3839 125.902 7.6318 183.169 -18.2011 20 0
66.8919 56.8502 131.803 7.48117 192.691 -17.9938 21 0
70.6081 56.3147 141.639 7.23013 202.205 -17.8278 22 0
74.3243 55.2473 149.508 7.02929 209.823 -17.6619 23 0
78.0405 54.7118 159.344 6.77824 219.353 -17.4133 24 0
81.7568 54.7082 169.18 6.42678 226.97 -17.2474 25 0
85.473 54.7046 177.049 6.27615 234.596 -17.0402 26 0
89.1892 55.7648 184.918 6.0251 244.11 -16.8742 27 0
92.9054 56.2931 192.787 5.82427 255.529 -16.6668 28 0
98.4797 57.3516 202.623 5.67364 265.059 -16.4181 29 0
100.338 57.8817 214.426 5.57322 276.477 -16.2106 30 0
107.77 58.4064 224.262 5.4728 289.793 -16.003 31 0
109.628 58.4046 238.033 5.37238 299.307 -15.837 32 0
115.203 58.9311 247.869 5.27197 306.924 -15.6712 33 0
118.919 58.9275 259.672 5.17155 316.423 -15.5878 34 0
122.635 58.392 267.541 5.02092 325.922 -15.5044 35 0
130.068 57.321 275.41 4.9205 331.627 -15.4214 36 0
135.642 57.3156 283.279 4.82008 339.229 -15.3381 37 0
143.074 56.2446 293.115 4.66946 344.927 -15.2964 38 0
148.649 55.7073 300.984 4.46862 352.529 -15.2132 39 0
156.081 55.7001 306.885 4.21757 358.218 -15.2128 40 0
161.655 55.6947 312.787 4.01674 365.82 -15.1295 41 0
169.088 55.6875 316.721 3.91632 371.51 -15.1291 42 0
178.378 56.2105 322.623 3.66527 379.112 -15.0459 43 0
185.811 56.7352 328.525 3.46444 388.603 -15.0039 44 0
191.385 57.2617 334.426 3.16318 396.197 -14.962 45 0
198.818 57.7864 338.361 2.91213 403.783 -14.9614 46 0
206.25 58.8431 342.295 2.7113 413.274 -14.9194 47 0
213.682 59.3678 346.23 2.46025 426.55 -14.9184 48 0
222.973 59.3588 354.098 2.20921 447.412 -14.9168 49 0
226.689 59.8871 356.066 2.05858 466.386 -14.8741 50 0
234.122 59.88 363.934 1.85774 483.455 -14.8728 51 0
243.412 60.4029 369.836 1.60669 504.318 -14.8712 52 0
247.128 59.8674 373.77 1.45607 519.49 -14.8701 53 0
250.845 59.3319 381.639 1.20502 536.559 -14.8688 54 0
256.419 59.3265 387.541 1.00418 549.835 -14.8678 55 0
261.993 58.7892 391.475 0.953975 56 0
269.426 57.7182 399.344 0.803347 57 0
271.284 57.1844 407.213 0.65272 58 0
276.858 56.1152 419.016 0.502092 59 0
282.432 55.5779 430.82 0.401674 60 0
288.007 53.9768 434.754 0.351464 61 0
293.581 52.9076 442.623 0.301255 62 0
297.297 51.3082 452.459 0.251046 63 0
304.73 49.7053 462.295 0.200837 64 0
308.446 48.106 474.098 0.200837 65 0
315.878 46.503 485.902 0.200837 66 0
319.595 45.4356 493.77 0.200837 67 0
321.453 43.8381 503.607 0.200837 68 0
325.169 42.7706 521.311 0.150628 69 0
330.743 41.7014 535.082 0.150628 70 0
334.459 40.634 546.885 0.150628 71 0
341.892 39.563 550.82 0.150628 72 0
349.324 39.0239 73 0
353.041 37.9564 74 0
358.615 36.3553 75 0
364.189 35.2861 76 0
371.622 33.6832 77 0
373.48 33.1494 78 0
377.196 32.6139 79 0
382.77 30.4809 80 0
390.203 28.8779 81 0
395.777 27.8087 82 0
403.209 25.6739 83 0
410.642 24.6029 84 0
416.216 23.0017 85 0
421.791 21.9325 86 0
429.223 20.3296 87 0
431.081 18.732 88 22.9
438.514 16.5972 89 22.9
445.946 15.5262 90 82.1
453.378 14.4551 91 82.1
458.953 12.854 92 90
464.527 11.7848 93 90
468.243 11.2493 94 90
475.676 9.64635 95 90
483.108 8.57533 96 180
488.682 7.50611 97 180
497.973 6.43329 98 180
503.547 5.89599 99 180
512.838 4.82317 100 180
518.412 4.28587 101 180
523.986 3.21665 102 180
531.419 3.20946 103 180
538.851 3.20227 104 180
542.568 2.13485 105 180
550 1.06383 106 180
107 180
108 180
109 180
110 180
111 180
112 180
113 180
114 180
115 180
116 180
117 180
118 180
119 180
120 180
121 180
122 180
123 180
124 180
125 180
126 180
127 180
128 180
129 180
130 180
131 180
132 180
133 180
134 180
135 180
136 180
137 180
138 0
139 0
140 0
141 0
142 0
143 0
144 0
145 0
146 0
147 0
148 0
149 0
150 0
151 0
152 0
153 0
154 0
155 0
156 0
157 0
158 0
159 0
160 0
161 0
162 0
163 0
164 0
165 0
166 15.4
167 15.4
168 15.4
169 15.4
170 41.7
171 41.7
172 41.7
173 41.7
174 49.5
175 49.5
176 49.5
177 49.5
178 49.5
179 49.5
180 35.3
181 35.3
182 32.4
183 32.4
184 32.4
185 32.4
186 32.4
187 32.4
188 32.4
189 32.4
190 42.4
191 42.4
192 42.4
193 42.4
194 42.4
195 42.4
196 42.4
197 42.4
198 42.4
199 42.4
200 51.9
201 51.9
202 51.9
203 51.9
204 51.9
205 51.9
206 51.9
207 51.9
208 51.9
209 51.9
210 54.2
211 54.2
212 54.2
213 54.2
214 54.2
215 54.2
216 54.6
217 54.6
218 -55.3
219 -55.3
220 -55.3
221 -55.3
222 -55.3
223 -55.3
224 -55.3
225 -55.3
226 -55.3
227 -55.3
228 -55.3
229 -55.3
230 -68.3
231 -68.3
232 -68.3
233 -68.3
234 -68.3
235 -68.3
236 -68.3
237 -68.3
238 -68.3
239 -68.3
240 -74.2
241 -74.2
242 -74.2
243 -74.2
244 -74.2
245 -74.2
246 -74.2
247 -74.2
248 -74.2
249 -74.2
250 -79.3
251 -79.3
252 -79.3
253 -79.3
254 -79.3
255 -79.3
256 -79.3
257 -79.3
258 -79.3
259 -79.3
260 -81.1
261 -81.1
262 -81.1
263 -81.1
264 -81.1
265 -81.1
266 -81.1
267 -81.1
268 -81.1
269 -81.1
270 -82.3
271 -82.3
272 -82.3
273 -82.3
274 -82.3
275 -82.3
276 -82.3
277 -82.3
278 -82.3
279 -82.3
280 -81.9
281 -81.9
282 -81.9
283 -81.9
284 -81.9
285 -81.9
286 -81.9
287 -81.9
288 -81.9
289 -81.9
290 -82.3
291 -82.3
292 -82.3
293 -82.3
294 -82.3
295 -82.3
296 -82.3
297 -82.3
298 -82.3
299 -82.3
300 -77.2
301 -77.2
302 -77.2
303 -77.2
304 -77.2
305 -77.2
306 -77.2
307 -77.2
308 -77.2
309 -77.2
310 -73.5
311 -73.5
312 -73.5
313 -73.5
314 -73.5
315 -73.5
316 -73.5
317 -73.5
318 -73.5
319 -73.5
320 -70.4
321 -70.4
322 -43.7
323 -43.7
324 -43.7
325 -43.7
326 -43.7
327 -43.7
328 -43.7
329 -43.7
330 -32.1
331 -32.1
332 -32.1
333 -32.1
334 -32.1
335 -32.1
336 -32.1
337 -32.1
338 -58.8
339 -58.8
340 62.8
341 62.8
342 63
343 63
344 63
345 63
346 63
347 63
348 63
349 63
350 91.9
351 91.9
352 94
353 94
354 94
355 94
356 94
357 94
358 94
359 94
360 78.5
361 78.5
362 78.5
363 78.5
364 78.5
365 78.5
366 78.5
367 78.5
368 78.5
369 78.5
370 72.5
371 72.5
372 72.5
373 72.5
374 72.5
375 72.5
376 72.5
377 72.5
378 72.5
379 72.5
380 64.1
381 64.1
382 57
383 57
384 -48.9
385 -48.9
386 -48.9
387 -48.9
388 -48.9
389 -48.9
390 -65.5
391 -65.5
392 -65.5
393 -65.5
394 -65.5
395 -65.5
396 -81.1
397 -81.1
398 -81.1
399 -81.1
400 -84.7
401 -84.7
402 -84.7
403 -84.7
404 -84.7
405 -84.7
406 -84.7
407 -84.7
408 -84.7
409 -84.7
410 -83.3
411 -83.3
412 -83.3
413 -83.3
414 -83.3
415 -83.3
416 -83.3
417 -83.3
418 -86
419 -86
420 -86
421 -86
422 -86
423 -86
424 -86
425 -86
426 -86
427 -86
428 93.5
429 93.5
430 102.6
431 102.6
432 103.4
433 103.4
434 105
435 105
436 88.8
437 88.8
438 88.8
439 88.8
440 88.8
441 88.8
442 88.8
443 88.8
444 88.8
445 88.8
446 88.8
447 88.8
448 88.8
449 88.8
450 88.8
451 88.8
452 88.8
453 88.8
454 88.8
455 88.8
456 88.8
457 88.8
458 88.8
459 88.8
460 88.8
461 88.8
462 88.8
463 88.8
464 88.8
465 88.8
466 88.8
467 88.8
468 88.8
469 88.8
470 88.8
471 88.8
472 88.8
473 88.8
474 88.8
475 88.8
476 88.8
477 88.8
478 88.8
479 88.8
480 88.8
481 88.8
482 88.8
483 88.8
484 88.8
485 88.8
486 88.8
487 88.8
488 88.8
489 88.8
490 88.8
491 88.8
492 88.8
493 88.8
494 88.8
495 88.8
496 88.8
497 88.8
498 0
499 0
500 0
501 0
502 0
503 0
504 0
505 0
506 0
507 0
508 0
509 0
510 0
511 0
512 0
513 0
514 0
515 0
516 0
517 0
518 0
519 0
520 0
521 0
522 0
523 0
524 0
525 0
526 0
527 0
528 0
529 0
530 0
531 0
532 0
533 0
534 0
535 0
536 0
537 0
538 0
539 0
540 0
541 0
542 0
543 0
544 0
545 0
546 0
547 0
548 0
549 0
550 0
Upvotes: 3
Views: 71
Answers (1)
John Hedengren
Reputation: 14346
Here are a few suggestions to help with troubleshooting and verifying the solution:
- Plug in the optimal trajectory with
STATUS=0to fix the MVs and see if the objective function decreases. If the solver reports a successful solution then the KKT Conditions are satisfied and the solution is optimal if the problem is convex. - If the problem is non-convex then try a multi-start method. A multi-start method initializes with different initial conditions to see if a different or better local solution is obtained.
- Try increasing the number of time steps or the accuracy of those time steps with
m.options.NODES=3. By default, Gekko uses 2 Nodes per time step. Here is more information on Nodes in APMonitor and Gekko. This grid-independence study is important to verify that there is no solution dependency on the time discretization.
These are some of the main things to try but there are other methods to verify and troubleshoot applications. Let us know if they help.
Upvotes: 1
Related Questions
- Change ODE initial condition to meet final target
- Optimal Control Problem with nonlinear equations
- How to use Gekko to solve for optimal control for a reusable reentry vehicle
- How can I force initial conditions for optimal control problem in Gekko?
- GEKKO optimal control with control variables
- Infrequent measurements in Gekko with extra simulation points
- Optimal control with time varying parameters and fixed control in Gekko
- Defining control variables and objective function in GEKKO
- Use Gekko and Python to fit a numerical ODE solution to data
- Simulating a vector ODE using python GEKKO