Unable to find explicit solution in Lagrangian optimization

I am trying to find the analytical solution to the following problem:

I tried solving it by coding the Lagrangian by hand and use solve, but Matlab prints the warning: "Unable to find explicit solution".

I used the following code:

syms e1 e2 p1 p2 rho gamma lambda
syms E H(e1,e2)

H(e1,e2) = (e1^rho +e2^rho)^(1/rho)

L(e1, e2, lambda) = p1*e1 +p2*e2 + lambda*(H(e1,e2)-E)
L_e1 = diff(L,e1) == 0
L_e2 = diff(L,e2) == 0
L_lambda = diff(L,lambda) == 0

system = [L_e1,L_e2,L_lambda]
[e1_s,e2_s,lambda_s]=solve(system,[e1 e2 lambda])

Do you know what I could do to solve this? Or is there a different and better way to find an analytical solution?

[Matt J]

Note that the problem can always be rewritten in the simpler form,

where x=e/E and P=p*E. This is assuming E is a known positive constant.

Answer 1

[Catalytic]

An analytical solution for 0<rho<1 is -

A=[1 0;
   0 1;
  -1 0;
  0 -1]*E;

[fval,i]=min(A*[p1;p2]);

e1=A(i,1);
e2=A(i,2);

[Catalytic]

You can see this graphically by plotting the constrained region. The region always has extreme points at (), so that's where the optimum must lie.

E=1;

for rho=[0.1:0.2:0.9]
    
    fimplicit(@(e1,e2)  abs(e1).^rho + abs(e2).^rho - E.^rho, [-1.5,1.5]); hold on
    
end

[Matt J]

I like it. And, in fact, because the extreme points lie at points where H(e1,e2) is not differentiable, it shows that you will never find the true solution with Lagrange multiplier analysis.

Answer 2

[Matt J]

If you make rho explicit, it seems to be able to find solutions. I doubt there would be a closed-form solution for general rho.

rho=2;
syms e1 e2 p1 p2 gamma lambda
syms E H(e1,e2)

H(e1,e2) = (e1^rho +e2^rho)

H(e1, e2) = 

L(e1, e2, lambda) = p1*e1 +p2*e2 + lambda*(H(e1,e2)-E^rho)

L(e1, e2, lambda) = 

L_e1 = diff(L,e1) == 0

L_e1(e1, e2, lambda) = 

L_e2 = diff(L,e2) == 0

L_e2(e1, e2, lambda) = 

L_lambda = diff(L,lambda) == 0

L_lambda(e1, e2, lambda) = 

system = [L_e1,L_e2,L_lambda]

system(e1, e2, lambda) = 

[e1_s,e2_s,lambda_s]=solve(system,[e1 e2 lambda])

e1_s = 

e2_s = 

lambda_s = 

[Fabian]

Thank you Matt! I think the problem arises for , is there a way to do the calculations for only?

[Torsten]

is there a way to do the calculations for rho in [0 1] only?

No. Even if rho is a numerical value, you won't get an explicit solution for e1 in most cases because you had to solve the below equation df=0 explicitly for e1:

syms e1 e2 rho E p1 p2
f = p1*e1 + p2*(E^rho-e1^rho)^(1/rho)

f = 

df = diff(f,e1)

df = 

solve(df==0,e1)

Warning: Unable to find explicit solution. For options, see help.

 
ans =
 
Empty sym: 0-by-1
 

[Matt J]

Even when it can be explicitly solved, the result isn't nice:

rho=sym(1/4);

syms e1 e2 p1 p2 gamma lambda
syms H(e1,e2)

H(e1,e2) = (e1^rho +e2^rho);

L(e1, e2, lambda) = p1*e1 +p2*e2 + lambda*(H(e1,e2)-1);
L_e1 = diff(L,e1) == 0;
L_e2 = diff(L,e2) == 0;
L_lambda = diff(L,lambda) == 0;

system = [L_e1,L_e2,L_lambda];
[e1_s,e2_s,lambda_s]=solve(system,[e1 e2 lambda])

Warning: Possibly spurious solutions.

e1_s = 

e2_s = 

lambda_s = 

[Walter Roberson]

You can eliminate the root() constructs, but the result is confusing.

rho=sym(1/4);

syms e1 e2 p1 p2 gamma lambda
syms H(e1,e2)

H(e1,e2) = (e1^rho +e2^rho);

L(e1, e2, lambda) = p1*e1 +p2*e2 + lambda*(H(e1,e2)-1);
L_e1 = diff(L,e1) == 0;
L_e2 = diff(L,e2) == 0;
L_lambda = diff(L,lambda) == 0;

system = [L_e1,L_e2,L_lambda];
[e1_s,e2_s,lambda_s]=solve(system,[e1 e2 lambda], 'maxdegree', 3)

Warning: Possibly spurious solutions.

e1_s = 

e2_s = 

lambda_s =