Good morning friends and Professor Hedengren, I am new to Python and even more to Gekko, first of all, I would like to know if my code in Gekko is correct, that is, according to what I physically want, considering that my equations are correct.
My model tries to optimize the variables M2 and l_M2 (or the combination of these 2 variables), in order to minimize in module the amplitude of movement (positive or negative) of my variable q1, my model receives inputs from the placed text file here, the model solution must respect the following:
- With the initial values of
M2andl_M2, the model is solved and the maximum amplitude (positive or negative) ofq1is obtained; - The input values do not vary over the horizon;
- In each iteration the value of the variable
c_m2must be updated according to the value ofM2andl_M2, and it must remain constant throughout the horizon.
In order to minimize the variable q1 I proposed two types of objectives, which I do not use simultaneously:
- Minimize
1000*q1**2; - Minimize the integral of
x1 = integral (0.5 q1 ** 2)dt evaluated from 0 to t, for which create an auxiliary variablex1.
Doubts to solve
- When solving the model, I realized that the value of
c_m2(at the initial point) is 0. Which is not correct, as it should be the same value as the following, so there is an error in my code , which I don't know. How to solve; - On the other hand, I would like to be able to compare the responses of the model with the initial values of the variables versus the response with the optimized values (as shown in the figure), but I cannot understand how to save my response with the initial values. Optimization check figure
- Is it correct to use
m.options.IMODE = 6in this case?
this is my code:
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()
