# Solve struggles to solve equation with variable constraint

I am trying to find {x,y} set satisfying first order conditions with the constraint:

0<x<=y<1


Here is my code, which does not give anything.

Dbuy = 1 - (x/y) - x*Log[y/x];
Dwait = (x/y) - x;
revx = D[(Dbuy*x + Dwait*CCC*y), x]
revy = D[(Dbuy*x + Dwait*CCC*y), y]
revsln=Solve[{revx==0, revy==0},{x,y}]


When I use FindRoot function instead of Solve for a fixed value of CCC (between 0 and 1, for example 0.4), it gives me numerical values (0.322, 0.581) if I add

{{x,0.01}, {y,0.01}}


But I want some form of closed-form solution for

0<=CCC<=1


What do you guys think?

• Well, do you have any reason to expect a closed-form solution exists at all? – AccidentalFourierTransform Oct 21 '18 at 20:16
• I feel that the root for the constraint 0<=x<=y<=1 has a closed-form. Is there a way to check it? Or is there a way to add a constraint to "Solve" function? Or is there a way to find a closed-form approximation? – Ovunc Oct 21 '18 at 20:20
• In 11.3 I get a "This system cannot be solved with the methods available to Solve." message. Probably means the system cannot be solved with the methods available to Solve, no? Even if I drop the floating CCC I get Root objects parametrized near to a value. – b3m2a1 Oct 21 '18 at 20:36
• Thanks b3m2a1, any suggestions? – Ovunc Oct 21 '18 at 23:09

Until someone comes up with a method that yields an analytical solution, you can use NArgMax:

argmax[c_?NumericQ] := NArgMax[{Dbuy*x + Dwait*c*y,
0 <= c <= 1 , 0 <= x <= 1, 0 <= y <= 1, x <= y}, {x, y}]

ParametricPlot[Evaluate @ argmax[c], {c, 0, 1},
PlotStyle -> Arrowheads[{.05, .05, .05, .05}],
PlotRangeClipping -> False,
PlotRange -> {{0, 1}, {0, 1}}, Frame -> True,
FrameLabel -> {{Style["y(c)", 16], None}, {Style["x(c)", 16], None}},
GridLines -> (List /@ argmax[1/2]), MeshFunctions -> {#3 &},
Mesh -> {{{0, Directive[PointSize[Large], Red]},
{1/2,  Directive[PointSize[Large], Orange]} ,
{1, Directive[PointSize[Large], Green]}}},
PlotLegends -> PointLegend[{Red, Orange, Green}, {"c = 0" , "c = 1/2", "c = 1"}],
Prolog -> {Gray, Arrowheads, Thin, Line[{{0, 0}, {1, 1}}]} ] /.
Line -> Arrow • are these the optimal (x,y) combinations for 0<=c<=1? – Ovunc Oct 21 '18 at 22:55
• @Ovunc, yes they are. – kglr Oct 21 '18 at 22:56
• This is very cool! Thank you. Is there any way to map at least a few $c$ values on this? – Ovunc Oct 21 '18 at 22:59
• @Ovunc, please see the updated version. – kglr Oct 21 '18 at 23:15
• Wow this is super cool! Thanks so much kglr! – Ovunc Oct 21 '18 at 23:16

Here's a way to get a sense for the kinds of things for want. We'll solve this (exactly!) for a large number of values for CCC, plot them, and extract meaning out of that.

Here's the basic code I'll fux with:

blurg[CCC_] :=
Block[
{
ccc = Rationalize[CCC],
Dwait,
revx,
revy,
x, y
},
Dbuy = 1 - (x/y) - x*Log[y/x];
Dwait = (x/y) - x;
revx = D[(Dbuy*x + Dwait*ccc*y), x];
revy = D[(Dbuy*x + Dwait*ccc*y), y];
Solve[{revx == 0, revy == 0}, {x, y}, Reals]
]

solns =
blurg /@ Rest@Subdivide[0, 1, 100];

N[{x, y} /. #, 50] & /@ solns //
ListPlot[#, PlotStyle -> Map[Hue, Range[0, 1, 1/Length@solns]]] & You can see you always have two solutions and they converge down to those end-points at CCC=1 and apparently diverge as CCC->0

You can also look at the symbolic solutions returned but they will be nasty. Lots of complicated Root objects

Here's another fun example:

solns =
blurg /@ Rest@Subdivide[0, .01, 25];

N[{x, y} /. #, 50] & /@ solns //
ListPlot[#, PlotStyle -> Map[Hue, Range[0, 1, 1/Length@solns]]] & This shows the divergence over that side of the solutions.

• Interesting work. I ran your code and saw the optimal (x,y) sets, but only numbers. I also saw the two-solution issue, but it is not a problem because x>y in the second solution of each c-value. – Ovunc Oct 21 '18 at 23:05
• I wonder if you can tell me a bit how this code works and how I can see the symbolic solutions. – Ovunc Oct 21 '18 at 23:06
• @Ovunc all of the solns give you symbolic solutions. I just provided a series of exact values for CCC and told Solve to solve it over the real numbers (that might not be necessary but I didn't really care to change it). I then extracted them with ReplaceAll and numericized them with N. – b3m2a1 Oct 21 '18 at 23:08
• What change should I make in the code to see the symbolic solutions? Are there more than two of them? – Ovunc Oct 21 '18 at 23:14

A straightforward solution might be found using NMinimize with constraint 0 < x <= y < 1:

sol[CCC_?NumericQ] := {x, y} /.
NMinimize[{1, {1 + CCC + x - (2 x)/y - CCC y - 2 x Log[y/x] == 0, -((x (x (-1 + y) + CCC y^2))/y^2) == 0,
0 < x <= y < 1}}  (* constraint *)
, {x, y}][]


Plot gives (after a while...)

Show[{ParametricPlot[ sol[c], {c, 0, 1}, PlotStyle -> Red,PlotRange -> {{0, 1}, {0, 1}}, AxesLabel -> {x[c], y[c]}],
RegionPlot[0 < x <= y < 1, {x, 0, 1}, {y, 0, 1}]}] the solution {x[c],y[c]} , 0<c<1 (same solution @kglr gave)

Although it may not be possible to obtain a symbolic solution for {x, y} as a function of CCC, it is straightforward to obtain a symbolic solution for {x, y, CCC} as a function of x/y. The derivation proceeds as follows:

s = Solve[{revx == 0, revy == 0} /. y -> a x, {x, CCC}] // Simplify
(* {{x -> (4 - a + Sqrt[8 - 8 a + a^2 + 8 Log[a]])/(4 a - 4 a Log[a]),
CCC -> -((-2 + a + Sqrt[8 - 8 a + a^2 + 8 Log[a]])/(2 a))},
{x -> (-4 + a + Sqrt[8 - 8 a + a^2 + 8 Log[a]])/(4 a (-1 + Log[a])),
CCC -> (2 - a + Sqrt[8 - 8 a + a^2 + 8 Log[a]])/(2 a)}} *)


To plot this result, it is convenient to obtain limits on a corresponding to those on CCC.

alim0 = Flatten[FindRoot[# == 0, {a, 1}] & /@ (CCC /. s)] // Values
(* {1., 3.51286} *)

alim1 = Chop[Flatten[Quiet@FindRoot[# == 1, {a, .8}] & /@ (CCC /. s)], 10^-7]
// Values
(* {0.692201, 1.} *)


The plot shown in earlier answers can be obtained by

Show[Plot[Null, {x, .25, 1}, PlotRange -> {Full, {.49, 1}}, ImageSize -> Large,
AxesLabel -> {x, y}, LabelStyle -> {Bold, Black, 12}],
ParametricPlot[{x, a x} /. First[s], {a, First@alim1, 1}],
ParametricPlot[{x, a x} /. Last[s], {a, 1, Last@alim0}]] And, the CCC dependence is shown by.

Show[Plot3D[Null, {x, .25, 1}, {y, .5, 1}, PlotRange -> {Full, Full, {0, 1}},
ImageSize -> Large, AxesLabel -> {x, y, CCC}, LabelStyle -> {Bold, Black, 12}],
ParametricPlot3D[{x, a x, CCC} /. First[s], {a, First@alim1, 1}],
ParametricPlot3D[{x, a x, CCC} /. Last[s], {a, 1, Last@alim0}]] Imposing the additional constraint given in the question, 0 < x <= y < 1 is achieved by using only the second solution, Last[s], which appears at the left of the two plots. The explicit dependence of {x, y} is, then,

ParametricPlot[Evaluate[{{CCC, x}, {CCC, a x}} /. Last[s]], {a, 1, Last@alim0},
ImageSize -> Large, AxesLabel -> {CCC, "x, y"},
LabelStyle -> {Bold, Black, 12}, AspectRatio -> 1/GoldenRatio] In response to a comment by the OP, below, the maximum value of the objective function can be determined and marked on a plot of the objective function (here with CCC == 0.1) can be accomplished by

dta = (FindMaximum[{Dbuy*x + Dwait*0.1*y, x <= y, 0 < x < 1, 0 < y < 1}, {x, y}]
// Reverse // Flatten) /. Rule[_, z_] -> z
(* {0.286834, 0.785105, 0.105362} *)

Show[Plot3D[Dbuy*x + Dwait*0.1*y, {x, 0, 1}, {y, 0, 1},
RegionFunction -> Function[{x, y}, x - y <= 0],
PlotRange -> {0, Automatic}, Mesh -> None, ImageSize -> Large,
AxesLabel -> {x, y, f}, LabelStyle -> {Black, Bold, 12}],
ListPointPlot3D[{dta}, PlotStyle -> {Red, Medium}]] • Interesting, it feels like this graph is the right of two. – Ovunc Oct 22 '18 at 12:14
• Can we show that the actual function is concave somehow? – Ovunc Oct 22 '18 at 12:31
• To which function or functions are you referring? There are several in your question and my answer. – bbgodfrey Oct 22 '18 at 12:41
• The objective function, (Dbuy*x + Dwait * ccc * y), is not concave. – Ovunc Oct 22 '18 at 13:05
• And thanks so much! I am working on your code right now! – Ovunc Oct 22 '18 at 13:08