How to find a proper algorithm to solve this optimal control problem?

25 views (last 30 days)
Ehsan Ranjbari
Ehsan Ranjbari on 13 Jan 2022
Commented: Ehsan Ranjbari on 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 Comments
Torsten
Torsten on 13 Jan 2022
Sorry, but I have no experience with numerical optimal control.
So I can't give you advise in this respect.

Sign in to comment.

Accepted Answer

Torsten
Torsten on 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 Comments
Ehsan Ranjbari
Ehsan Ranjbari on 21 Jan 2022
I managed to reformulate the problem and actually simplified it. Now the problem reduces to this: https://www.mathworks.com/matlabcentral/answers/1630720-solve-optimal-control-problem-with-free-final-time-using-matlab?s_tid=srchtitle
I would love to hear your suggestion on solving this.
Thanks in Advance

Sign in to comment.

More Answers (1)

Walter Roberson
Walter Roberson on 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
f = 
Dfu = diff(f,u)
Dfu = 
string(Dfu)
ans = "((7*x)/1000 + 3/500)*(70*(2/(5*x))^(1/10) - 100)*(x - 1) - (2*u)/25 - 7/100"
solu = simplify(solve(Dfu, u))
solu = 
Dfx = diff(f,x)
Dfx = 
string(Dfx)
ans = "(70*(2/(5*x))^(1/10) - 100)*((3*u)/500 + x*((7*u)/1000 + 3/1000) + 1/200) + ((7*u)/1000 + 3/1000)*(70*(2/(5*x))^(1/10) - 100)*(x - 1) - (14*(x - 1)*((3*u)/500 + x*((7*u)/1000 + 3/1000) + 1/200))/(5*x^2*(2/(5*x))^(9/10))"
  3 Comments
Ehsan Ranjbari
Ehsan Ranjbari on 14 Jan 2022
Thank you both for the comments.
Yes. The problem is exactly the on @Torsten mentioned above.
This problem is generally an economic one and the aim is to maximize the that integral (market share).
I want to share my effort on using Hamiltonian:
%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)
syms f ;
f = ((p - c0*((x0/x)^z))*Dx) - (a + (b*u) + (c*u^2));
%% Hamiltonian %lambda_0= 1 (Normal case)
syms H p1 ;
H = f + p1*Dx;
%% Costate equations
Dp1 = -diff(H,x);
%% solve for control u
du = diff(H,u);
sol_u = solve(du,u);
%% Substitute u to state equation
Dx = subs(Dx, u, sol_u);
f = subs(f, u, sol_u);
%% Convert symbolic objects to strings for using 'dsolve'
eq1 = strcat('Dx=',char(Dx));
eq2 = strcat('Dp1=',char(Dp1));
sol_h = dsolve(eq1,eq2);
%% Use boundary conditions to determine the coefficients
% conA1 = 'x(0) = 0.4';
% conA2 = 'x(31) = 1';
% sol_a = dsolve(eq2,conA1,conA2);
after running this I get the warning for eq1 which is:
% Warning: Unable to find symbolic solution.
I think it is better to use the numerical solutions and see what happens.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!