Maximized value of a function as an additional constraint on a variable

Question

I want to maximize a function $f(b, h)$ with respect to its arguments and subject to some additional constraints. If these constraints are satisfied, $f$ is increasing in $h$. I then want to find an upper bound $u$ on $h$ such that the maximized value of $f$ remains smaller or equal than some threshold level.

More precisely, I have the following function:

f[h_, b_] := .1 (1 + b) + (0.3 b (-b + h))/(-1 + h) 

A standard numerical maximization of $f$, NMaximize[{f[h, b], 1 <= b <= h <= u}, {h, b}] yields {1.50833, {h -> 10., b -> 6.5}} for $u = 10$.

What I would like to have is another function which maximizes the value of the upper bound $u$ subject to a constraint on the maximized function value.

For example:

If the maximized function value of $f$ was to remain below 1.4, $u$ would have to be less than 9.18125, since NMaximize[{f[h, b], 1 <= b <= h <= 9.18125}, {h, b}] yields {1.4, {h -> 9.18125, b -> 5.95418}}.

The following function $f_m(u)$ captures the first part of what I want to do:

fm[u_?NumericQ] := Block[{h, b}, {h, b} /. NMaximize[{f[h, b], 1 <= b <= h <= u}, {h, b}]] 

Calling $f_m(u)$ for different $u$ generates different values of the maximized function that are increasing in $u$, as can be seen from tabulating the results by Table[fm[u][[2]], {u, 2, 10}]:

{0.508, {h -> 2., b -> 1.167}}
{0.604, {h -> 3., b -> 1.833}}
{0.725, {h -> 4., b -> 2.500}}
{0.852, {h -> 5., b -> 3.167}}
{0.981, {h -> 6., b -> 3.834}}
{1.112, {h -> 7., b -> 4.500}}
{1.244, {h -> 8., b -> 5.169}}
{1.376, {h -> 9., b -> 5.835}}
{1.508, {h -> 10., b -> 6.500}}

What I would now like to do is something like this (which does not work, however): NMaximize[{fm[u], fm[[u]][[2]][[1]] < 1.4}, {u}]. That is, I want to select the highest value of $u$ such that the value of $f_m(u)$ does not exceed the threshold of 1.4.

What also does not work is NMaximize[{f[h, b], 1 <= b <= h <= u, f[h, b] < 1.4}, {h, b}] yielding {1.4, {h -> 9.5519, b -> 5.01912}} because it picks one point on the level curve at which $f$ equals 1.4. I am looking for a specific point on that curve, however, namely the point where, given $b$ is maximized, $h$ approaches the upper bound $u$.

I have tried several variants of what has been suggested here: Combined numerical minimization and maximization, adapted to two successive maximizations for $b$ and $h$. I did not get this to work, however, with the constraints on $b$ included in the inner NMaximize command.

I also tried to work around the bi-variate maximization by using using the envelope theorem to express the maximum of $b$ in terms of $h$ and then compute the threshold. This also did not work, however, because it disregards the constraints on $b$ and $h$, but $f$ approaches infinity as $h$ goes to 1 for $b \le h$ , so NMaximize finds the maximum at $h = 1$.

Any suggestions are greatly appreciated!

Answer 1

You can use FindRoot!

First, take a look at the contour plot of $f$:

ContourPlot[f[h, b], {h, 1, 6}, {b, 0, 6},
 PlotPoints -> 100,
 Contours -> {0.1, 0.3, 0.5, 0.7, 0.9},
 ContourLabels -> True,
 ContourShading -> None,
 ContourStyle -> GrayLevel@0.7,
 FrameLabel -> {h, b}]

Note that there is a saddle point at $(h, b) = (1.75, 1)$. The location can be determined using Solve:

Solve[f[h, b] == 0.5, b]

{{b -> 1.}, {b -> 1.33333 (-1. + h)}}

The intersection of these two solutions is the saddle point:

Solve[%[[1, 1, 2]] == %[[2, 1, 2]], h]

{{h -> 1.75}}

If $b$ is restricted to be greater than or equal to 1, then for any $h$ to the left of the saddle point, $f_\textrm{max}(h, b) = 0.5$, while for $h$ to the right, $f_\textrm{max}(h, b)$ is greater than 0.5 and increases monotonically with $h$, following the line of steepest ascent evident on the contour plot.

One can find the maximum value of $f$ for a given $u$ as you have indicated, but including the extra Method -> "SimulatedAnnealing" seems to help Mathematica deal with the kink at the saddle point:

fm[u_?NumericQ] := NMaxValue[
  {f[h, b], 1 <= b <= h <= u}, {h, b}, 
  Method -> "SimulatedAnnealing"]

This is a perfectly good numeric function, so you can throw it into FindRoot:

um[fmax_?NumericQ] := u /. FindRoot[fm[u] == fmax, {u, 10}]

This will work for $f_\textrm{max} > 0.5$ (there is no solution below 0.5). If you need solutions very close to 0.5, you may need to increase the precision all around. At $f_\textrm{max} = 0.5$, as there are a range of $u$-values for which $f_m = 0.5$, it is best to give the function the value we know it should take, namely the $h$-coordinate of the saddle point:

um[0.5] = h /. First@Solve[Equal @@ Solve[f[h, b] == 0.5, b][[All, 1, 2]], h];

It's plottable:

Plot[
 um[fmax], {fmax, 0.5, 10},
 Axes -> False,
 Frame -> True,
 FrameLabel -> {Subscript[f, max], u}
]