Btukfyl

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;