我的问题在于约束状态变量（寻找 23、24、25）。我目前正在使用 ODE45 来求解方程，并使用 fmincon 来找到最佳控制变量。你将如何解决这个问题？到目前为止，这是我的代码：

function ster_6

optionsopt = optimset('Display','iter-detailed','Algorithm','sqp', 'TolFun', 10^(-8), 'TolX', 10^(-8), 'MaxFunEvals', 10^6, 'MaxIter', 50);

u_lb = [ 0.0];

u_ub = [ 50.0]

u0 = (u_lb+u_ub)/2;

[rozw, dokladnosc] = fmincon(@model_procesu,u0,[],[],[],[],u_lb,u_ub,[],optionsopt)

end

function J = model_procesu(u)

tspan_1 = [0.0  140.0];
x0_1 = [ 1.5    0.0    0.0    7.0 ];


options_ode = odeset('RelTol',1e-6,'AbsTol',1e-6);

[tsol_1,xsol_1] = ode45(@(t,x) penicillin_production_problem( t,x,u(1) ), tspan_1, x0_1, options_ode);

J = -(xsol_1(end,2) * xsol_1(end,4));


plot(tsol_1,xsol_1)
title('Trajektorie zmiennych stanu')

ylabel('Stężenie')

xlabel('Czas')
legend('x1_1(t)','x2_1(t)','x3_1(t)','x4_1(t)')

end

function dx = penicillin_production_problem(t,x,u)

dx = zeros(4,1);


h1 = 0.11 * (x(3)/(0.006*x(1)+x(3)));

h2 = 0.0055 * (x(3)/(0.0001 + x(3) * (1 + (10 * x(3)))));

dx(1) = h1 * x(1) - u * ( x(1) / (500*x(4)) );

dx(2) = h2 * x(1) - 0.01*x(2) - u * ( x(2) / (500*x(4)));

dx(3) = (-((h1 * x(1)) / 0.47)) - (h2 * (x(1) / 1.2)) - (x(1) * ((0.029*x(3)) / (0.0001 + x(3)))) + ((u/x(4)) * (1-(x(3)/500)));

dx(4) = u/500;

end


```