在 GAMS 中使用辛普森积分规则的问题
Problem using Simpson's integration rule in GAMS
我使用 GAMS 编写了一个简单的代码，它使用梯形积分确定滑翔机的最大范围。我想用 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),
*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;
这是使用辛普森的有缺陷的：
$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;
我做的就是照着书做的
使用非线性规划进行最优控制和估计的实用方法 第 141 页和第 142 页之间。由于我的控制未知，y_hat 只是以下总和的平均值y_k+1 和 y_k，然后，我在这些点定义了变量 D 和 L，然后计算了 y_k+1 - y_k 它是如何在第 141 页中被引用的。
但是，现在我看到的不是第一个代码中显示的变量，而是某种奇怪的循环。 This is my propper answer with trapezoid rule
and this is my defective solution with Simpson's method.
非常欢迎所有关于我的错误或错误所在的建议。
感谢阅读。
经过一段时间的尝试，我发现是许可证问题导致了这些问题。对代码进行简单的更改即可使其正常工作。
$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;
我使用 GAMS 编写了一个简单的代码，它使用梯形积分确定滑翔机的最大范围。我想用 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),
*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;
这是使用辛普森的有缺陷的：
$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;
我做的就是照着书做的
使用非线性规划进行最优控制和估计的实用方法 第 141 页和第 142 页之间。由于我的控制未知，y_hat 只是以下总和的平均值y_k+1 和 y_k，然后，我在这些点定义了变量 D 和 L，然后计算了 y_k+1 - y_k 它是如何在第 141 页中被引用的。
但是，现在我看到的不是第一个代码中显示的变量，而是某种奇怪的循环。 This is my propper answer with trapezoid rule and this is my defective solution with Simpson's method.
非常欢迎所有关于我的错误或错误所在的建议。 感谢阅读。
经过一段时间的尝试，我发现是许可证问题导致了这些问题。对代码进行简单的更改即可使其正常工作。
$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;