問題是通過改變kst,x1,x5和xo，在(-8e-4到2e-4)的范圍內找到x3的最佳(最大)值。
x5=5 %Input 2（輸入2是一個狀態變數，在執行時可能在4到15范圍內變化
優化時可能在4到15的范圍內變化)
kst=1 %輸入3（輸入3是速率常數，它可以在0.1到2之間變化）。
xo=4 %Input 4 (Input 4是一個狀態變數，可以在4到10的范圍內變化)。
x1=1e-7 %輸入1可以在1e-9到1e-6之間變化。
- 腳本檔案
function rest = Scrpt1（t，X）
x2 = X(1)。
x3 = X(2)。
%Parameters[/span
if t<15
x1 = 1e-7; %輸入1可以從1e-9到1e-6變化
else x1 = 0;
end
x5=5 %Input 2 (Input 2是一個狀態變數，在進行優化時可以在4到15的范圍內變化)
kst=1 %輸入3（輸入3是速率常數，它可以在0.1到2的范圍內變化）。
xo=4 %輸入4（輸入4是一個狀態變數，可以在4到10的范圍內變化）。
k1 = 6e7;
km1 = 0.20;
km4 = 0.003;
k3 = 2500.00;
k4 = km4/9;
km3 = km1;
LAP=1.5; LAP=1.5
%的微分方程。
dx2dt = km1*x3 km3*LAP - k1*x1*x2 km4*x3 - k4*x2;
dx3dt = k1*x1*x2 - km1*(x3 x5 xo) - k3*x3*kst;
rest = [dx2dt; dx3dt]。
end。
- 用于ODE求解的函式檔案 。
options = odeset('InitialStep',0. 0001,'RelTol',1e-09）。)
[T,Y]=ode15s(@Scrpt1,[0 60],[9e-13,0], options);
X3= Y(:,2)。
plot(T,X3)
如何使用fmincon或任何其他優化求解器來解決所提到的尋找x3的最大值的優化問題。對于x5,kst,xo,x1的哪個值，我們可以得到最大的x3?
uj5u.com熱心網友回復：
首先你必須添加你想優化的值作為你的耦合衍射方程的引數：
function rest = Scrpt1（t。 X,X_opt）
x5=X_opt(1)。
kst=X_opt(2)。
xo=X_opt(3)。
x1=X_opt(4)。
x2 = X(1);
x3 = X(2)。
%Parameters[/span
if t>=15
x1 =0。
結束。
k1 = 6e7;
km1 = 0.20;
km4 = 0.003;
k3 = 2500.00;
k4 = km4/9;
km3 = km1;
LAP=1.5;
%的微分方程。
dx2dt = km1*x3 km3*LAP - k1*x1*x2 km4*x3 - k4*x2;
dx3dt = k1*x1*x2 - km1*(x3 x5 xo) - k3*x3*kst;
rest = [dx2dt; dx3dt]。
end。
然后你必須撰寫一個你想最小化的包裝函式。因為你想使x3最大化，所以你必須在目標值上加一個減號。
函式 max_X3=fun（X_opt）
tspan=[0 60] 。
y0=[9e-13,0] 。
options = odeset('InitialStep',0.0001,'RelTol',1e-09) 。
[~,y] = ode15s(@(t,y) Scrpt1(t,y,X_opt), tspan, y0,options)。
max_X3=-max（y(:,2）。
結束。
最后你可以像這樣使用fmincon：
% x5, kst, xo, x1
initial_search_point=[5, 1, 4, 1-7]
lower_bounds=[4, 0.1, 4, 1-9]
upper_bounds=[15, 2, 10, 1e-6]
fmincon(@fun,initial_search_point,[],[],[], lower_bounds,upper_bounds)
uj5u.com熱心網友回復：
下面是一個解決方案，在
import numpy as np
from gekko import GEKKO
n = 121; t = np.linspace(0,60, n)
m = GEKKO(remote=False)
m.time = t
k1 = 6e7; km1 = 0.20; km4 = 0.003。
k3 = 2500.00; k4 = km4/9;
km3 = km1; LAP=1.5; k4 = km4/9; km3 = km1; LAP=1.5
x5 = m.FV(value=5,lb=4,ub=15); x5.STATUS = 1.
kst = m.FV(value=1,lb=0.1,ub=2); kst.STATUS =1.
xo = m.FV(value=4,lb=4,ub=10); xo.STATUS = 1.
x1 = m. FV(value=[1e-17 if t[i]< 15 else 0 for i in range（n）】。
lb=1-9,ub=1-6)
x2,x3 = m.Array(m.Var,2)
x3.值 = -0.00032
x3.lower = -8e-4[/span]。
x3.upper = 2e-4[/span]。
m.Equations([x2.dt()==km1*x3 km3*LAP-k1*x1*x2 km4*x3-k4*x2,
x3.dt()==k1*x1*x2-km1*(x3 x5 xo)-k3*x3*kst])
m.Maximize(x3)
m.options.SOLVER =1
m.options.IMODE = 6
m.solve()
import matplotlib.pyplot as plt
plt.plot(m.time,x3)
plt.show()
x2和x3的初始條件沒有定義。
狀態變數的數量。 483
總方程的數量。- 480 總方程數：480。
松弛變數的數量: - 0
---------------------------------------
自由度：3
----------------------------------------------
用APOPT求解器進行動態控制
----------------------------------------------
迭代目標收斂性
0 4.99217E 02 2.99935E-011 6.07645E-02 4.31439E-05
2 3.25294E-02 3.04712E-05
3 3.41027E-02 8.96081E-054 3.31615E-02 2.48287E-06
5 3.31615E-02 2.22045E-16
63.31615E-022.22045E-16
成功的解決方案
---------------------------------------------------
解算器 : APOPT (v1.0)
求解時間 : 2.189999245629E-002 sec
目標：3.316154172805905E-002
成功的解決方案
---------------------------------------------------
轉載請註明出處，本文鏈接：https://www.uj5u.com/gongcheng/309713.html
標籤：