Problem with Solve (Lagrange multiplier)

Question

I'm trying to solve a simple maximization problem, which I can easily solve by hand using the Lagrangian method.

I want to maximize UH subject to BC=0.

 UH[CH_, IH_] := CH^α (γH V1 + γH IH)^(1 - α)
BC[CH_, IH_] := CH + IH + ρ τ (γH V1 + γH IH) - YH

$Assumptions = CH > 0 && CH ∈ Reals && IH > 0 && IH ∈ Reals && V1 > 0 && V1 ∈ Reals && YH > 0 && YH ∈ Reals && α > 0 && α < 1 && α ∈ Reals && γH > 
1 && γH ∈ Reals && lambda ∈ Reals ;

Other threads recommend using Maximize, but

Maximize[{UH[CH, IH], BC[CH, IH] == 0}, {CH, IH}]

does not return any result.

I have also tried to implement the Lagrangian with

L = UH[CH, IH] - lambda (BC[CH, IH])
Solve[{{D[L, CH] == 0, D[L, IH] == 0, D[L, lambda] == 0}, CH, IH}, lambda]

which returns the error message that this is not a quantified system of equations and inequalities.

I am probably overlooking something simple (I'm new to Mathematica), but I have now spent two hours on this simple problem and would thus really appreciate any help.

Answer 1

The Lagrange multiplier method is just a way of expressing the requirement that the gradients of the target function and the constraint function are linearly dependent. Therefore, you can implement the method directly by requiring that the determinant of the Jacobian should vanish. Here I do this:

Solve[
 BC[CH, IH] == 0 && 
  Simplify[Det[D[{UH[CH, IH], BC[CH, IH]}, {{CH, IH}}]] == 0], {CH, IH}]

(*
==> {{CH -> V1 α + YH α, 
  IH -> -((-YH + V1 α + YH α + 
     V1 γH ρ τ)/(1 + γH ρ τ))}}
*)

The D operation acts on the vector of both functions, and therefore gives the Jacobian matrix. I combine this with the constraint equation and then solve for the two unknowns in one step. Simplify is inserted before doing Solve to apply the $Assumptions as early as possible to eliminate undesirable solutions.

The difference to the approach using Maximize is that this method doesn't actually verify whether it's a maximum. But that's precisely why the above code works even without specifying numerical values of the parameters.

Answer 2

The constraints of the optimization problem are both equality constraints (BC) and inequality constraints (CH>0 and IH>0). The Lagrangian for the maximization problem is defined as:

Lg[CH_, IH_, x_, y_, z_] := UH[CH, IH] - x BC[CH, IH] + y CH + z IH

the multiplier x is associated with the equality constraint while the multipliers y and z are associated with the non-negativity constraints on the variables vars={CH,IH}.

A full characterization of the solution should in principle consider the first order conditions focs = Thread[D[Lg[CH, IH, x, y, z], {{CH, IH, x}}] == 0] along with the complementary slackness constraints cscs = Thread[{y, z} vars == 0], the original non-negativity constraints oineqcs = Thread[vars > 0] and the non-negativity of the multipliers of the inequality constraints nnoineqcs = Thread[{y, z} >= 0].

Now, consider the fact that for the complementary slackness conditions to hold, the corresponding multipliers must equal 0 because of the strict non-negativity of the variables (oineqcs). Because of this, the critical points of the system can be obtained by solving the system of first order conditions (the derivatives of the Lagrangian wrt the variables vars and the equality constraint multiplier x):

Solve[
  (* enforce the complementary slackness requirement and solve *)
  focs /. Thread[{y, z} -> 0]//Simplify, Join[vars, {x}]
  (* simplify using assumptions about parameters *)
 ] // Simplify[#, Assumptions -> gineqs] &

{{CH -> (V1 + YH) α, 
    IH -> (YH - YH α - V1 (α + γH ρ τ))/(1 + γH ρ τ),
    x -> -(-1 + α) γH (α/(γH - α γH))^α (1/(1 + γH ρ τ))^(1 - α)}, 
   {CH -> Undefined, IH -> Undefined, x -> Undefined}, 
   {CH -> Undefined, IH -> Undefined, x -> Undefined}}

The (first) solution above is a maximizer if the determinant of the bordered hessian of the Lagrangian is positive ie

(V1 + YH)^(-1 + α) α^(-1 + α) γH (-(((V1 + YH) (-1 + α) γH)/(1 + γH ρ τ)))^-α (1 + γH ρ τ) > 0

Note: The expression above is derived as Det[bh] where bh is the bordered hessian (derived below) and evaluated at the maximizer

 bh = With[{h = D[Lg[CH, IH, x, y, z], {{CH, IH}, 2}], c = D[BC[CH, IH], {{CH, IH}}]},
   Transpose[Prepend[Transpose[Prepend[h, c]], Prepend[c, 0]]]]

(also, gineqs = {V1 > 0, YH > 0, α > 0, α < 1, γH > 1} are the inequalities associated with the various parameters)