How to find a proper algorithm to solve this optimal control problem?
107 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Ehsan Ranjbari
am 13 Jan. 2022
Kommentiert: Ehsan Ranjbari
am 21 Jan. 2022
Hi everyone!
I am trying to find a way to solve this optimal control problem in MATLAB. The function is too complex and the using Hamiltonian in MATLAB couldn't help.The problem describes as below:
p = 100;
a = 0;
b = 0.07;
c = 0.04;
r = 0.005;
z = 0.1;
c0 = 70;
x0 = 0.4;
alpha = 0.005;
beta = 0.006;
gamma = 0.003;
delta = 0.007;
Dx = (alpha + beta*u + (gamma + delta*u)*x)*(1-x); % State Equation
f = ((p - c0*((x0/x)^z))*Dx) - (a + (b*u) + (c*u^2)); % Function inside the integral (Cost function)
% x(t0) = 0.4, x(tf) = free, t0 = 0. tf = 31
Note that the aim is to maximize the function f.
I tried to use fmincon and still the function is too complex to get an answer.
Thanks!
3 Kommentare
Torsten
am 13 Jan. 2022
Sorry, but I have no experience with numerical optimal control.
So I can't give you advise in this respect.
Akzeptierte Antwort
Torsten
am 14 Jan. 2022
This should give you a start:
%Optimal advertising expenditure in monopoly
%% Constants
p = 100;
a = 0;
b = 0.07;
c = 0.04;
r = 0.005;
z = 0.1;
c0 = 70;
x0 = 0.4;
alpha = 0.005;
beta = 0.006;
gamma = 0.003;
delta = 0.007;
%% State equation (g)
syms x u p1
Dx = (alpha + beta*u + (gamma + delta*u)*x)*(1-x);
%% Cost function inside the integral (f)
f = ((p - c0*((x0/x)^z))*Dx) - (a + (b*u) + (c*u^2));
%% Hamiltonian %lambda_0= 1 (Normal case)
H = f + p1*Dx;
%% Costate equations
Dp1 = -diff(H,x);
%% solve for control u
du = diff(H,u);
sol_u = solve(du,u);
f = subs(f,u,sol_u)
Dp1 = subs(Dp1,u,sol_u)
rhs = [f;Dp1];
% Turn to numerical computation
fun = matlabFunction(rhs)
tmesh = linspace(0,31,150);
guess = @(x)[0.4*(1-x/31)+x/31;1]
solinit = bvpinit(tmesh,guess);
bvpfcn = @(t,y)fun(y(2),y(1));
bcfcn = @(ya,yb)[ya(1)-0.4;yb(1)-1];
sol = bvp4c(bvpfcn, bcfcn, solinit)
2 Kommentare
Torsten
am 14 Jan. 2022
Since you want to maximize, you'll have to take
f = -(((p - c0*((x0/x)^z))*Dx) - (a + (b*u) + (c*u^2)));
I guess.
Weitere Antworten (1)
Walter Roberson
am 13 Jan. 2022
You did not say what you wanted to optimzie with respect to. If you wanted to optimize with respect to u, then see solu below.
If you wanted to optimize with respect to x (in terms of u) then I will need to do more testing.
syms x u
p = 100;
a = 0;
b = 0.07;
c = 0.04;
r = 0.005;
z = 0.1;
c0 = 70;
x0 = 0.4;
alpha = 0.005;
beta = 0.006;
gamma = 0.003;
delta = 0.007;
Dx = (alpha + beta*u + (gamma + delta*u)*x)*(1-x); % State Equation
f = ((p - c0*((x0/x)^z))*Dx) - (a + (b*u) + (c*u^2)); % Function inside the integral (Cost function)
f
Dfu = diff(f,u)
string(Dfu)
solu = simplify(solve(Dfu, u))
Dfx = diff(f,x)
string(Dfx)
3 Kommentare
Torsten
am 13 Jan. 2022
Bearbeitet: Torsten
am 13 Jan. 2022
I think the problem is
Find u(t) such that
integral_{t=0}^{t=tf} ((p - c0*((x0/x(t))^z))*dx/dt) - (a + (b*u(t)) + (c*u(t)^2)) dt
is maximized under the constraint
dx/dt = (alpha + beta*u(t) + (gamma + delta*u(t))*x(t))*(1-x(t))
x(0)=0, x(tf) = free
Siehe auch
Kategorien
Mehr zu Conversion Between Symbolic and Numeric finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!