我使用 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;
我所做的就是照着书
Practical Methods for Optimal Control and Estimation Using Nonlinear Programming在第 141 和 142 页之间。由于我的控制是未知的，因此 y_hat 只是 y_k+1 和 y_k 之和的平均值，因此，我在这些点定义了变量 D 和 L，并且然后计算 y_k+1 - y_k 在第 141 页中的建议方式。
但是，现在我看到的不是第一个代码中显示的变量，而是某种奇怪的循环。这是我使用梯形规则的正确答案， 这是我使用辛普森方法的有缺陷的解决方案。
所有关于我的错误或错误的地方的建议都非常受欢迎。谢谢阅读。