Добрый день!
Прежде всего, я некоторое время искал похожие вопросы, но (вероятно, из-за моей неопытности) я не нашел ничего похожего на то, что я собираюсь спросить.
Я впервые использую matlab для решения подобных проблем, поэтому я не уверен, что делать.Краткое объяснение:
Я делаю проект для своего курса «Оптимальное управление»: мне нужно воспроизвести результаты статьи о трудоустройстве, и я застрял на графиках.У меня есть следующие данные:
- пять переменных функций (U (t), T (t), R (t), V1 (t) и V2 (t))
- четыре функции управления (u1 (t), u2 (t), u3 (t), u4 (t))
- ограничения на управляющие переменные (каждое u должно быть в диапазоне от 0 до 1)
- начальные значения для U, T, R, V1 и V2 (при t = 0, в частности, V1 и V2 постоянны во времени)
- конечные значения для коэффициентов λ в гамильтониане
(примечание: для элементов управления я уже нашел оптимальное выражение, которое находится в такой форме: ui = min {1, max {0, "expression"}}. При необходимости я могу также указать четыревыражения, которые немного пренебрегли синтезом)
По совету профессора я попытался использовать fmincon, который теоретически должен дать мне информацию, которая мне нужна, чтобы построить какой-то результат, используя только функцию стоимости задачи.Но в этом случае у меня есть некоторые проблемы, связанные со временем в расчетах.Ниже код, который я использовал для fmincon:
syms u
%note: u(5) corresponds to U(t), but this is the only way I've found to get
%a result, the other u(i) are in ascending order (u(1) = u1 and so on...)
g = @(u) 30*u(5) + (20/2)*(u(1))^2 + (20/2)*(u(2))^2 + (10/2)*(u(3))^2 + (40/2)*(u(4))^2;
%initial guesses
u0 = [0 0 0 0 100000]; %
A = [];
b = [];
Aeq = [];
beq = [];
lb = 0.0 * ones(1,2,3,4);
ub = 1.0 * ones(1,2,3,4);
[x,fval,output,lambda] = fmincon(g, u0, A, b, Aeq, beq, lb, ub);
С этим кодом я получаю (очевидно) только одно значение для каждой переменной в качестве результата, и, поскольку я не нашел ни одного метода, связанного со временем,как я уже говорил, я начинаю искать другие стратегии решения.
Я обнаружил, что ode45 - это средство для решения дифференциальных уравнений, в которое уже включена «итерация времени», поэтому я попытался написать код, чтобы получитьэто работает с моей проблемой.Я взял все уравнения из бумаги и поместил их в вектор, как показано в примерах математических работ, и это мой файл matlab:
syms u1(t) u2(t) u3(t) u4(t)
syms U(t) T(t) R(t) V1(t) V2(t)
syms lambda_u lambda_t lambda_r lambda_v1 lambda_v2
%all the parameters provided by the paper
delta = 500;
alpha1 = 0.004;
alpha2 = 0.005;
alpha3 = 0.006;
gamma1 = 0.001;
gamma2 = 0.002;
phi1 = 0.22;
phi2 = 0.20;
delta1 = 0.09;
delta2 = 0.05;
k1 = 0.000003;
k2 = 0.000002;
k3 = 0.0000045;
%these two variable are set constant
V1 = 200;
V2 = 100;
%weight values for the cost function (only A1 is used in this case, but I left them all since the unused ones are irrelevant)
A1 = 30;
A2 = 20;
A3 = 20;
A4 = 10;
A5 = 40;
%ordering the unknowns in an array
x = [U T R u1 u2 u3 u4];
%initial conditions, ordered as the x vector (for the ui are guesses)
y0 = [100000 2000 1000 0 0 0 0];
%system set up
f = @(t,x) [delta - (1 + x(4))*k1*x(1)*V1 - (1 + x(5))*k2*x(1)*V2 - alpha1*x(1) + gamma1*x(2) + gamma2*x(3);...
(1 + x(4))*k1*x(1)*V1 - k3*x(2)*V2 - alpha2*x(2) - gamma1*x(2);...
(1 + x(5))*k2*x(1)*V2 - alpha3*x(3) - gamma2*x(3) + k3*x(2)*V2;...
alpha2*x(2) + gamma1*x(2) + (1 + x(6))*phi1*x(1) + k3*x(2)*V2 - delta1*V1;...
alpha3*x(3) + gamma2*x(3) + (1 + x(7))*phi2*x(1) - delta2*V2;...
-A1 + (1 + x(4))*k1*V1*(lambda_u - lambda_t) + (1 + x(5))*k2*V2*(lambda_u - lambda_r) + lambda_u*alpha1 - lambda_v1*(1 + x(6))*phi1 - lambda_v2*(1 + x(7))*phi2;...
-lambda_u*gamma1 + (alpha2 + gamma1)*(lambda_t - lambda_v1) + k3*V2*(lambda_t - lambda_r - lambda_v1);...
-lambda_u*gamma2 + (alpha3 + gamma2)*(lambda_r - lambda_v2);...
(1 + x(4))*k1*x(1)*(lambda_u - lambda_t) + lambda_v1*delta1;...
(1 + x(5))*k2*x(1)*(lambda_u -lambda_r) + k3*x(2)*(lambda_t - lambda_r - lambda_v1) + lambda_v2*delta2];
%using ode45 to solve over the chosen time interval
[t,xa] = ode45(f,[0 10],y0);
С этим кодом я получаю следующую ошибку:
Error using odearguments (line 95)
@(T,X)[DELTA-(1+X(4))*K1*X(1)*V1-(1+X(5))*K2*X(1)*V2-ALPHA1*X(1)+GAMMA1*X(2)+GAMMA2*X(3);(1+X(4))*K1*X(1)*V1-K3*X(2)*V2-ALPHA2*X(2)-GAMMA1*X(2);(1+X(5))*K2*X(1)*V2-ALPHA3*X(3)-GAMMA2*X(3)+K3*X(2)*V2;ALPHA2*X(2)+GAMMA1*X(2)+(1+X(6))*PHI1*X(1)+K3*X(2)*V2-DELTA1*V1;ALPHA3*X(3)+GAMMA2*X(3)+(1+X(7))*PHI2*X(1)-DELTA2*V2;-A1+(1+X(4))*K1*V1*(LAMBDA_U-LAMBDA_T)+(1+X(5))*K2*V2*(LAMBDA_U-LAMBDA_R)+LAMBDA_U*ALPHA1-LAMBDA_V1*(1+X(6))*PHI1-LAMBDA_V2*(1+X(7))*PHI2;-LAMBDA_U*GAMMA1+(ALPHA2+GAMMA1)*(LAMBDA_T-LAMBDA_V1)+K3*V2*(LAMBDA_T-LAMBDA_R-LAMBDA_V1);-LAMBDA_U*GAMMA2+(ALPHA3+GAMMA2)*(LAMBDA_R-LAMBDA_V2);(1+X(4))*K1*X(1)*(LAMBDA_U-LAMBDA_T)+LAMBDA_V1*DELTA1;(1+X(5))*K2*X(1)*(LAMBDA_U-LAMBDA_R)+K3*X(2)*(LAMBDA_T-LAMBDA_R-LAMBDA_V1)+LAMBDA_V2*DELTA2]
returns a vector of length 10, but the length of initial conditions vector is 7. The vector returned by
@(T,X)[DELTA-(1+X(4))*K1*X(1)*V1-(1+X(5))*K2*X(1)*V2-ALPHA1*X(1)+GAMMA1*X(2)+GAMMA2*X(3);(1+X(4))*K1*X(1)*V1-K3*X(2)*V2-ALPHA2*X(2)-GAMMA1*X(2);(1+X(5))*K2*X(1)*V2-ALPHA3*X(3)-GAMMA2*X(3)+K3*X(2)*V2;ALPHA2*X(2)+GAMMA1*X(2)+(1+X(6))*PHI1*X(1)+K3*X(2)*V2-DELTA1*V1;ALPHA3*X(3)+GAMMA2*X(3)+(1+X(7))*PHI2*X(1)-DELTA2*V2;-A1+(1+X(4))*K1*V1*(LAMBDA_U-LAMBDA_T)+(1+X(5))*K2*V2*(LAMBDA_U-LAMBDA_R)+LAMBDA_U*ALPHA1-LAMBDA_V1*(1+X(6))*PHI1-LAMBDA_V2*(1+X(7))*PHI2;-LAMBDA_U*GAMMA1+(ALPHA2+GAMMA1)*(LAMBDA_T-LAMBDA_V1)+K3*V2*(LAMBDA_T-LAMBDA_R-LAMBDA_V1);-LAMBDA_U*GAMMA2+(ALPHA3+GAMMA2)*(LAMBDA_R-LAMBDA_V2);(1+X(4))*K1*X(1)*(LAMBDA_U-LAMBDA_T)+LAMBDA_V1*DELTA1;(1+X(5))*K2*X(1)*(LAMBDA_U-LAMBDA_R)+K3*X(2)*(LAMBDA_T-LAMBDA_R-LAMBDA_V1)+LAMBDA_V2*DELTA2]
and the initial conditions vector must have the same number of elements.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in test (line 62)
[t,xa] = ode45(f,[0 10],y0);
Для которого я не могу найти решение, так как я использовал все начальные значения, приведенные в статье.Единственные значения, которые я оставил, являются окончательными значениями лямбда-коэффициентов, так как они являются окончательными значениями, и я не уверен, что их можно использовать.В этом случае я также не могу понять, где я должен поставить границы для контрольной переменной.
Для полноты я также предоставлю ссылку на рассматриваемый документ: https://www.ripublication.com/ijss17/ijssv12n3_13.pdf
Можете ли вы помочь мне понять, что я могу сделать, чтобы решить свои проблемы?
PS: я знаю, что это довольно плохой код, но я основываюсь на базовых уроках по математике;конечно, это должно быть реорганизовано и упорядочено в различных файлах (один для функции стоимости и один для ограничений, например), но сначала я хотел бы понять, где проблема, и затем я приведу все в симпатичную форму.
Спасибо большое!