Я сделал простой код, используя GAMS, который определяет максимальную досягаемость планера, используя интеграцию трапеции.Я хочу воссоздать ту же программу с интеграцией Симпсона, однако не могу понять результаты.
Это функциональный код с правилом трапеции:
$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_это просто среднее из суммы y_k + 1 и y_k, затем я определил переменные D и L в этих точках, а затем вычислил y_k + 1 - y_k, как это предлагается на странице 141.
Однако вместо того, чтобы видеть переменные, отображаемые как в первом коде, теперь я вижу какой-то странный цикл. Это мой правильный ответ с правилом трапеции , и это моё ошибочное решение по методу Симпсона .
Все рекомендации о том, где мои ошибки или ошибки очень приветствуются.Спасибо за чтение.