I have made a simple code using GAMS which determines the maximum reach of a glider using trapeziod integration. I want to recreate the same program with SImpson's integration, however, I cannot understand the results.
This is the functional code with the trapezoid rule:
$set n 50
set j /0*%n%/;
sets
jlast(j)
jnotlast(j);
jlast(j)$(ord(j)=card(j))=yes;
jnotlast(j)=not jlast(j);
scalar
n number of intervals /%n%/
m mass /5000/
S surface /21.55/
CD0 drag /0.023/
k ni idea /0.073/
hmax initial height /1000/
g gravity /9.81/
density density /1.225/
variable
gamma(j),
CL(j),
D(j),
CD(j),
L(j),
*x(j),
*y(j),
objective;
positive variable
x(j),
y(j),
v(j),
step;
equation
diffx(j),
diffy(j),
valueD(j),
valueL(j),
obj;
diffx[j]$(jnotlast(j)).. x[j+1]-x[j] =e=0.5*step*(v(j+1)*cos(gamma(j+1)) + v(j)*cos(gamma(j)) );
diffy[j]$(jnotlast(j)).. y[j+1]-y[j] =e=0.5*step*(v(j+1)*sin(gamma(j+1)) + v(j)*sin(gamma(j)) );
valueD[j].. m*g*sin(gamma(j))=e=0.5*density*S*v(j)*v(j)*(CD0+k*CL(j)*CL(j));
valueL[j].. m*g*cos(gamma(j))=e=0.5*density*S*v(j)*v(j)*CL(j);
obj .. objective =e= x('%n%');
x.fx('0') = 1.0e-12;
y.fx('0') = 1000;
y.fx('%n%') = 1.0e-12;
CL.up(j) =1.4;
y.up (j) = 1000;
gamma.up(j) = pi*0.5;
v.lo(j) = 1.0e-12;
y.lo(j) = 1.0e-12;
CL.lo(j) = 0;
gamma.lo(j) = 0;
model brahstron1 /all/;
option
nlp=ipopt;
solve brahstron1 using nlp maximize objective;
And this is the defective one using Simpson:
$set n 50
set j /0*%n%/;
sets
jlast(j)
jnotlast(j);
jlast(j)$(ord(j)=card(j))=yes;
jnotlast(j)=not jlast(j);
scalar
n number of intervals /%n%/
m mass /5000/
S surface /21.55/
CD0 drag /0.023/
k ni idea /0.073/
hmax initial height /1000/
g gravity /9.81/
density density /1.225/
variable
gamma(j),
CL(j),
D(j),
CD(j),
L(j),
gamma_med(j),
CL_med(j),
D_med(j),
CD_med(j),
L_med(j),
objective;
positive variable
x(j),
y(j),
v(j),
x_med(j),
y_med(j),
v_med(j),
step;
equation
diffx(j),
diffy(j),
diffx_central(j),
diffy_central(j),
valueD(j),
valueL(j),
valueD_central(j),
valueL_central(j),
obj;
diffx[j]$(jnotlast(j)).. x[j+1]-x[j] =e=(1/6)*step*(v(j+1)*cos(gamma(j+1)) + v(j)*cos(gamma(j)) + 4*v_med(j+1)*cos(gamma_med(j+1)) );
diffy[j]$(jnotlast(j)).. y[j+1]-y[j] =e=(1/6)*step*(v(j+1)*sin(gamma(j+1)) + v(j)*sin(gamma(j)) + 4*v_med(j+1)*sin(gamma_med(j+1)) );
diffx_central[j]$(jnotlast(j)).. x_med[j+1] =e=0.5*(x(j+1)+x(j));
diffy_central[j]$(jnotlast(j)).. y_med[j+1] =e=0.5*(y(j+1)+y(j));
valueD[j].. m*g*sin(gamma(j))=e=0.5*density*S*v(j)*v(j)*(CD0+k*CL(j)*CL(j));
valueL[j].. m*g*cos(gamma(j))=e=0.5*density*S*v(j)*v(j)*CL(j);
valueD_central[j].. m*g*sin(gamma_med(j))=e=0.5*density*S*v_med(j)*v_med(j)*(CD0+k*CL_med(j)*CL_med(j));
valueL_central[j].. m*g*cos(gamma_med(j))=e=0.5*density*S*v_med(j)*v_med(j)*CL_med(j);
obj .. objective =e= x('%n%');
x.fx('0') = 1.0e-12;
y.fx('0') = 1000;
y.fx('%n%') = 1.0e-12;
CL.up(j) =1.4;
CL_med.up(j) =1.4;
y.up (j) = 1000;
y_med.up (j) = 1000;
gamma.up(j) = pi*0.5;
gamma_med.up(j) = pi*0.5;
v.lo(j) = 1.0e-12;
v_med.lo(j) = 1.0e-12;
y.lo(j) = 1.0e-12;
y_med.lo(j) = 1.0e-12;
CL.lo(j) = 0;
CL_med.lo(j) =0;
gamma.lo(j) = 0;
gamma_med.lo(j) = 0;
model brahstron1 /all/;
* Invoke the LGO solver option for solving this nonlinear programming
option
nlp=ipopt;
solve brahstron1 using nlp maximize objective;
What I did was to follow the book
Practical Methods for Optimal Control and Estimation Using Nonlinear Programming between pages 141 and 142. Since my control is unknown the y_hat are simply the average of the sum of y_k+1 and y_k, then, I defined the variables D and L at these points and then calculated y_k+1 - y_k how it is sugested in page 141.
However, instead of seeing the variables displayed as in the first code, now I see some kind of weird loop. This is my propper answer with trapezoid rule and this is my defective solution with Simpson's method.
All recomendations on where my error or errors are are extremely welcome. Thanks for reading.
After triyng for some time, I have found out that it is a licence problem what is causing these issues. A simple change in the code enables it to work as it should.
$set n 10
set j /0*%n%/;
sets
jlast(j)
jnotlast(j);
jlast(j)$(ord(j)=card(j))=yes;
jnotlast(j)=not jlast(j);
scalar
n number of intervals /%n%/
m mass /5000/
S surface /21.55/
CD0 drag /0.023/
k ni idea /0.073/
hmax initial height /1000/
g gravity /9.81/
density density /1.225/
variable
gamma(j),
CL(j),
D(j),
CD(j),
L(j),
gamma_med(j),
CL_med(j),
D_med(j),
CD_med(j),
L_med(j),
objective;
positive variable
x(j),
y(j),
v(j),
x_med(j),
y_med(j),
v_med(j),
step;
equation
diffx(j),
diffy(j),
diffx_central(j),
diffy_central(j),
valueD(j),
valueL(j),
valueD_central(j),
valueL_central(j),
obj;
diffx[j]$(jnotlast(j)).. x[j+1]-x[j] =e=(1/6)*step*(v(j+1)*cos(gamma(j+1)) + v(j)*cos(gamma(j)) + 4*v_med(j+1)*cos(gamma_med(j+1)) );
diffy[j]$(jnotlast(j)).. y[j+1]-y[j] =e=(-1)* (1/6)*step*(v(j+1)*sin(gamma(j+1)) + v(j)*sin(gamma(j)) + 4*v_med(j+1)*sin(gamma_med(j+1)) );
diffx_central[j]$(jnotlast(j)).. x_med[j+1] =e=0.5*(x(j+1)+x(j)+(step/8)*(v_med(j)*cos(gamma_med(j)))-(v_med(j+1)*cos(gamma_med(j+1))));
diffy_central[j]$(jnotlast(j)).. y_med[j+1] =e=0.5*(y(j+1)+y(j)+(step/8)*(v_med(j)*sin(gamma_med(j)))-(v_med(j+1)*sin(gamma_med(j+1))));
valueD[j].. m*g*sin(gamma(j))=e=0.5*density*S*v(j)*v(j)*(CD0+k*CL(j)*CL(j));
valueL[j].. m*g*cos(gamma(j))=e=0.5*density*S*v(j)*v(j)*CL(j);
valueD_central[j].. m*g*sin(gamma_med(j))=e=0.5*density*S*v_med(j)*v_med(j)*(CD0+k*CL_med(j)*CL_med(j));
valueL_central[j].. m*g*cos(gamma_med(j))=e=0.5*density*S*v_med(j)*v_med(j)*CL_med(j);
obj .. objective =e= x('%n%');
x.fx('0') = 1.0e-12;
y.fx('0') = 1000;
y.fx('%n%') = 1.0e-12;
CL.up(j) =1.4;
CL_med.up(j) =1.4;
y.up (j) = 1000;
y_med.up (j) = 1000;
gamma.up(j) = pi*0.5;
gamma_med.up(j) = pi*0.5;
gamma.lo(j) = 0;
gamma_med.lo(j) = 0;
v.lo(j) = 1.0e-12;
v_med.lo(j) = 1.0e-12;
y.lo(j) = 1.0e-12;
y_med.lo(j) = 1.0e-12;
CL.lo(j) = 0;
CL_med.lo(j) =0;
gamma.lo(j) = 0;
gamma_med.lo(j) = 0;
model brahstron1 /all/;
* Invoke the LGO solver option for solving this nonlinear programming
option
nlp=ipopt;
solve brahstron1 using nlp maximize objective;