1
$\begingroup$

When I add a non-binding constraint to a maximization problem the result changes. I don't understand why.

Below you can find the code that I am using.

For the maximization A, H = 6837.66 For the maximization B (equal to A + constraint of H less or equal to 14000), H = 5716.08

If the solution in A is already below 14000, the solution should not change. Why is this happening?

Thank you!

$Assumptions -> {{β, θ, α, r, p, ϵ, σ, δ} > 0, Element [{H, β, θ, α, r, p, ϵ, σ, δ}, Reals], {β, θ, α, r, p, ϵ, σ, δ} < 1, H > 0, ϵ > σ, α < r};

Y = 1000;
β = 0.95;
θ = 10;
α = 0.05;
r = 0.05;
ϵ = 1.2;
σ = 0.2;
δ = 0.3;
p = 0.51;

g1[H_] := Max[0, ((((1 - δ)*Y)/(r + α)) - H)]*((((1 - δ)*Y)/(r + α)) - H)^(-1)

h1[H_] :=  Min[0, ((((1 - δ)*Y)/(r + α)) - H)]*((((1 - δ)*Y)/(r + α)) - H)^(-1)

A = Maximize[{(1 + β * (1 - 
          p))*(((Y - r*H)^(1 - 1/ϵ) + (θ*H)^(1 - 
             1/ϵ))^((1 - 1/σ)/(1 - 
            1/ϵ))) /(1 - 1/σ) +  β*p *
     g1[H] *(((Y - (r + α)*H)^(1 - 1/ϵ) + (θ*
             H)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 
            1/ϵ)) ) /(1 - 1/σ) + β*p *
     h1[H] *(((δ*Y)^(1 - 1/ϵ) + (θ*5000)^(1 - 
             1/ϵ))^((1 - 1/σ)/(1 - 
            1/ϵ)) ) /(1 - 1/σ), 
   H ⩾ 0, (Y - r*H) ⩾ 0}, {H}]



B = Maximize[{(1 + β * (1 - 
          p))*(((Y - r*H)^(1 - 1/ϵ) + (θ*H)^(1 - 
             1/ϵ))^((1 - 1/σ)/(1 - 
            1/ϵ))) /(1 - 1/σ) +  β*p *
     g1[H] *(((Y - (r + α)*H)^(1 - 1/ϵ) + (θ*
             H)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 
            1/ϵ)) ) /(1 - 1/σ) + β*p *
     h1[H] *(((δ*Y)^(1 - 1/ϵ) + (θ*5000)^(1 - 
             1/ϵ))^((1 - 1/σ)/(1 - 
            1/ϵ)) ) /(1 - 1/σ), 
   H ⩾ 0, (Y - r*H) ⩾ 0, 
   H <= 14000}, {H}]

```
Improve this question
$\endgroup$
4
  • 1
    $\begingroup$ In principle it shouldn't change, but in practice constrained vs. unconstrained optimization are two very different problems. The algorithms used behind the scenes to solve them may also be very different, maybe using different choices of starting values, and possibly leading to different results. How different are your results in the two cases? $\endgroup$
    – MarcoB
    Commented Apr 5, 2021 at 14:10
  • $\begingroup$ I need to simulate different solutions using different parameter values. When I create a vector for the solutions using option B the results stop making sense (as a vector). I know that the option A is the correct one, but as I need to compute thousands of solutions, I need to trust the function. At this moment, I don't feel comfortable at all with Maximize. $\endgroup$
    – Ana Sá
    Commented Apr 5, 2021 at 21:55
  • $\begingroup$ I wonder if it has to do with the discontinuity in your function. Your function has a division by zero when $H=\frac{(1-\delta ) Y}{\alpha +r}$. For your particular set of parameters that happens when $H=7000$. You might need to perform Maximize twice: once with $H<\frac{(1-\delta ) Y}{\alpha +r}$ and another with $H>\frac{(1-\delta ) Y}{\alpha +r}$. Then take the maximum of the two results. $\endgroup$
    – JimB
    Commented May 5, 2021 at 15:44
  • $\begingroup$ And setting $H$ to something just less than $\frac{Y-\delta Y}{\alpha +r}$ seems to be the maximum in this case. $\endgroup$
    – JimB
    Commented May 5, 2021 at 16:33

3 Answers 3

Reset to default
0
$\begingroup$

The option RandomSearch and the multiplication by 10^24 do the job.

Element[{H, \[Beta], \[Theta], \[Alpha], r, 
p, \[Epsilon], \[Sigma], \[Delta]},  Reals], {\[Beta], \[Theta], \[Alpha], r, 
p, \[Epsilon], \[Sigma], \[Delta]} < 1, H > 0, \[Epsilon] > \[Sigma], \[Alpha] < r};

Y = 1000;\[Beta] = 0.95;\[Theta] = 10;\[Alpha] = 0.05;r = 0.05;
\[Epsilon] = 1.2;\[Sigma] = 0.2;\[Delta] = 0.3;p = 0.51;
g1[H_] :=  Max[0, ((((1 - \[Delta])*Y)/(r + \[Alpha])) - 
 H)]*((((1 - \[Delta])*Y)/(r + \[Alpha])) - H)^(-1)
h1[H_] :=  Min[0, ((((1 - \[Delta])*Y)/(r + \[Alpha])) - 
 H)]*((((1 - \[Delta])*Y)/(r + \[Alpha])) - H)^(-1)

A = NMaximize[{10^24*((1 + \[Beta]*(1 - 
        p))*(((Y - r*H)^(1 - 1/\[Epsilon]) + (\[Theta]*H)^(1 - 
           1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 
          1/\[Epsilon])))/(1 - 1/\[Sigma]) + \[Beta]*p*
   g1[H]*(((Y - (r + \[Alpha])*H)^(1 - 1/\[Epsilon]) + (\[Theta]*
           H)^(1 - 1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 
          1/\[Epsilon])))/(1 - 1/\[Sigma]) + \[Beta]*p*
   h1[H]*(((\[Delta]*Y)^(1 - 1/\[Epsilon]) + (\[Theta]*5000)^(1 - 
           1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 
          1/\[Epsilon])))/(1 - 1/\[Sigma])), 
   H \[GreaterSlantEqual] 0, (Y - r*H) \[GreaterSlantEqual] 0}, {H}, 
  Method -> {"RandomSearch", "SearchPoints" -> 250}]

{-3.26868, {H -> 7000.}}

B = NMaximize[{10^24*((1 + \[Beta]*(1 - 
        p))*(((Y - r*H)^(1 - 1/\[Epsilon]) + (\[Theta]*H)^(1 - 
           1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 
          1/\[Epsilon])))/(1 - 1/\[Sigma]) + \[Beta]*p*
   g1[H]*(((Y - (r + \[Alpha])*H)^(1 - 1/\[Epsilon]) + (\[Theta]*
           H)^(1 - 1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 
          1/\[Epsilon])))/(1 - 1/\[Sigma]) + \[Beta]*p*
   h1[H]*(((\[Delta]*Y)^(1 - 1/\[Epsilon]) + (\[Theta]*5000)^(1 - 
           1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 
          1/\[Epsilon])))/(1 - 1/\[Sigma])), 
   H \[GreaterSlantEqual] 0, (Y - r*H) \[GreaterSlantEqual] 0, H <= 14000}, {H}, 
  Method -> {"RandomSearch", "SearchPoints" -> 250}]

{-3.26868, {H -> 7000.}}

Improve this answer
$\endgroup$
1
  • $\begingroup$ Taking "SearchPoints" -> 1000 in the above, we obtain close (but not identical) results without the multiplication by 10^24, $\endgroup$
    – user64494
    Commented May 7, 2021 at 12:15
0
$\begingroup$

Let us build the plot

Plot[(1 + \[Beta]*(1 - 
    p))*(((Y - r*H)^(1 - 1/\[Epsilon]) + (\[Theta]*H)^(1 - 
       1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 1/\[Epsilon])))/(1 - 
  1/\[Sigma]) + \[Beta]*p*   g1[H]*(((Y - (r + \[Alpha])*H)^(1 - 1/\[Epsilon]) + (\[Theta]*
       H)^(1 - 1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 
      1/\[Epsilon])))/(1 - 1/\[Sigma]) + \[Beta]*p* h1[H]*(((\[Delta]*Y)^(1 - 1/\[Epsilon]) + (\[Theta]*5000)^(1 - 
       1/\[Epsilon]))^((1 - 1/\[Sigma])/(1 - 1/\[Epsilon])))/(1 - 
  1/\[Sigma]), {H, 5000, 7000}]

enter image description here

I think the question is closed.

Addition. See the plot on {H,5000,15000} enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ The plot confirms that the maximization in A is the correct one. The problem is that I need to trust the maximize function and I need to create a vector for the solutions with different parameters and the constraint in B. $\endgroup$
    – Ana Sá
    Commented Apr 5, 2021 at 21:50
0
$\begingroup$

I think the issue is that your function is not continuous at $H= \frac{Y-\delta Y}{\alpha +r}$ as at that value of $H$ there is a division by zero.

If $H\leq \frac{Y-\delta Y}{\alpha +r}$, then g1[H]=1 and h1[H]=0. Otherwise, g1[H]=1 and h1[H]=0. That gives two functions:

f1 = ((1 + (1 - p) β) ((-H r + Y)^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/
  (1 - 1/σ) + (p β ((Y - H (r + α))^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/
  (1 - 1/σ)

f2 = (p β ((Y δ)^(1 - 1/ϵ) + 5000^(1 - 1/ϵ) θ^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/(1 - 1/σ) +
  ((1 + (1 - p) β) ((-H r + Y)^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/(1 - 1/σ)

Maximize both functions (with the appropriate restrictions on $H$), and then take the maximum result:

Y = 1000;
β = 0.95;
θ = 10;
α = 0.05;
r = 0.05;
ϵ = 1.2;
σ = 0.2;
δ = 0.3;
p = 0.51;
sol1 = Maximize[{f1, H <= (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0}, H]
(* {-3.26868*10^-24, {H -> 7000.}} *)
sol2 = Maximize[{f2, H > (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0}, H]
(* {-4.45483*10^-24, {H -> 14852.5}} *)
max = If[sol1[[1]] > sol2[[1]], sol1, sol2]
(* {-3.26868*10^-24, {H -> 7000.}} *)

This approach seems more stable when additional restrictions are added:

sol1 = Maximize[{f1, H <= (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0, H < 14000}, H]
(* {-3.26868*10^-24, {H -> 7000.}} *)
sol2 = Maximize[{f2, H > (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0, H < 14000}, H]
(* {-4.46379*10^-24, {H -> 14000.}} *)
max = If[sol1[[1]] > sol2[[1]], sol1, sol2]
(* {-3.26868*10^-24, {H -> 7000.}} *)

As a check on varying parameter values consider some sort of a Manipulate:

Manipulate[
 f1 = ((1 + (1 - p) β) ((-H r + Y)^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/(1 - 1/σ) + (p β ((Y - H (r + α))^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/(1 - 1/σ);
 f2 = (p β ((Y δ)^(1 - 1/ϵ) + 5000^(1 - 1/ϵ) θ^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/(1 - 1/σ) + ((1 + (1 - p) β) ((-H r + Y)^(1 - 1/ϵ) + (H θ)^(1 - 1/ϵ))^((1 - 1/σ)/(1 - 1/ϵ)))/(1 - 1/σ);
 sol1 = Maximize[{f1, H <= (Y - Y δ)/(r + α), H > 0, (Y - r*H) ⩾ 0}, H];
 sol2 = Maximize[{f2, H > (Y - Y δ)/(r + α),  H > 0, (Y - r*H) ⩾ 0}, H]; 
 max = If[sol1[[1]] > sol2[[1]], sol1, sol2];
 hmin = 0.8 H /. sol1[[2]];
 hmax = 1.2 H /. sol2[[2]];
 {ymin, ymax} = 
  MinMax[{f1 /. H -> hmin, f1 /. H -> (Y - Y δ)/(r + α),
     f2 /. H -> (Y - Y δ)/(r + α), f2 /. H -> hmax}];
 y0 = ymin - 0.95 (ymax - ymin);
 y1 = ymax + 1.05 (ymax - ymin);
 Show[Plot[f1, {H, hmin, (Y - Y δ)/(r + α)}, PlotRange -> {{hmin, hmax}, {y0, y1}}],
  Plot[f2, {H, (Y - Y δ)/(r + α), hmax}, PlotRange -> {{hmin, hmax}, {y0, y1}}],
  ListPlot[{{{H /. max[[2]], max[[1]]}}, {{H /. sol1[[2]], sol1[[1]]}},
    {{H /. sol2[[2]], sol2[[1]]}}}, PlotStyle -> {{PointSize[0.03], Red}, Blue, Blue}]],
 {{Y, 1000}, 1, 1000, Appearance -> "Labeled"},
 {{α, 0.05}, 0.02, 1, Appearance -> "Labeled"},
 {{r, 0.05}, 0.01, 1, Appearance -> "Labeled"},
 {{δ, 0.3}, 0.01, 1, Appearance -> "Labeled"},
 {{β, 0.95}, 0, 1, Appearance -> "Labeled"},
 {{θ, 10}, 0, 100, Appearance -> "Labeled"},
 {{ϵ, 1.2}, 0, 100, Appearance -> "Labeled"},
 {{σ, 0.2}, 0, 1, Appearance -> "Labeled"},
 {{p, 0.51}, 0, 1, Appearance -> "Labeled"},
 TrackedSymbols :> {Y, α, r, δ, β, θ, σ, p}]

Maxima for both parts of the function

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks for the accept. But when I bugged you, it was because you said you had thousands of these to do and because so much time had passed, I assumed you had given up on this forum. Someone else might still have a better and maybe more stable answer. $\endgroup$
    – JimB
    Commented May 7, 2021 at 3:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.