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

Question

Ehsan Ranjbari 2022-1-13

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1628245-how-to-find-a-proper-algorithm-to-solve-this-optimal-control-problem

评论： Ehsan Ranjbari 2022-1-21

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 个评论
显示 1更早的评论隐藏 1更早的评论

Ehsan Ranjbari 2022-1-13

Thanks for your kind reply.

I would also want to know if it is a good idea to first linearize the cost function and then using the Hamiltonian?

thanks in advance.

Torsten 2022-1-13

Sorry, but I have no experience with numerical optimal control.

So I can't give you advise in this respect.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-1-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1628245-how-to-find-a-proper-algorithm-to-solve-this-optimal-control-problem#answer_874595

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 个评论
显示无隐藏无

Torsten 2022-1-14

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.

Ehsan Ranjbari 2022-1-21

Dear @Torsten

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

请先登录，再进行评论。

Answer 2

Walter Roberson 2022-1-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1628245-how-to-find-a-proper-algorithm-to-solve-this-optimal-control-problem#answer_874185

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 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2022-1-13

编辑：Torsten 2022-1-13

@Walter Roberson

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

Ehsan Ranjbari 2022-1-14

@Walter Roberson and @Torsten.

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);