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?

  • 1
    Well, do you have any reason to expect a closed-form solution exists at all? – AccidentalFourierTransform Oct 21 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 at 20:20
  • 1
    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 at 20:36
  • Thanks b3m2a1, any suggestions? – Ovunc Oct 21 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[0], Thin, Line[{{0, 0}, {1, 1}}]} ] /. 
 Line -> Arrow

enter image description here

  • are these the optimal (x,y) combinations for 0<=c<=1? – Ovunc Oct 21 at 22:55
  • @Ovunc, yes they are. – kglr Oct 21 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 at 22:59
  • @Ovunc, please see the updated version. – kglr Oct 21 at 23:15
  • Wow this is super cool! Thanks so much kglr! – Ovunc Oct 21 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],
   Dbuy,
   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]]] &

enter image description here

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]]] &

enter image description here

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 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 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 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 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}][[2]]

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}]}]

enter image description here

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}]]

enter image description here

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}]]

enter image description here

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]

enter image description here

Addendum

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}]]

enter image description here

  • Interesting, it feels like this graph is the right of two. – Ovunc Oct 22 at 12:14
  • Can we show that the actual function is concave somehow? – Ovunc Oct 22 at 12:31
  • To which function or functions are you referring? There are several in your question and my answer. – bbgodfrey Oct 22 at 12:41
  • The objective function, (Dbuy*x + Dwait * ccc * y), is not concave. – Ovunc Oct 22 at 13:05
  • And thanks so much! I am working on your code right now! – Ovunc Oct 22 at 13:08

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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