1
$\begingroup$

I am having issues with the code:

CC = p0^2/(p0 - L)
Dbuy = (1 - CC)*(L/p0) + (CC - p0)*Log[CC/p0]
Dwait = (1 - p1)*(L/(p1 - p0))
Dlock = (1 - p1)*(1 - L/p0 - L/(p1 - p0)) + (p1 - CC)*(1 - L/p0) - L*Log[(p1 - p0)/(CC - p0)]
RevLocking1 = Table[NMaximize[{x*Dbuy*p0 + x*Dwait*0.4*p1 + x*Dlock*(p0 + 0.4*L), x*(Dbuy + Dwait + Dlock) <= 1,  0 <= L <= p0 <= p1 <= 1}, {p0, p1, L}, MaxIterations -> 10000], {x, 1, 6, 0.25}]

The code works, but it ignores the constraint:

0 <= L <= p0 <= p1 <= 1 

Since it reaches values like:

L -> 0.303712, p0 -> -0.0156877, p1 -> 0.900987

Any ideas?

$\endgroup$
3
  • $\begingroup$ Possible duplicate: mathematica.stackexchange.com/questions/59706/… $\endgroup$ – Michael E2 Oct 19 '18 at 21:21
  • $\begingroup$ I used the "Re" trick, but didn't really work. Do you have a better suggestion? $\endgroup$ – Ovunc Oct 19 '18 at 23:21
  • $\begingroup$ @ovunc: You should use && in between constraints. You are using ,, which is not correct. NMaximize[{obj, const1&&const2&&...},{vars}] $\endgroup$ – Tugrul Temel Oct 19 '18 at 23:49
1
$\begingroup$

It might be that the starting values are too far away from those that maximize the function. Rather than using Table one could put everything into a Do loop and use the previous estimates as the starting values using FindMaximum rather than NMinimize. Using (or stealing) some of the code from @TugrulTemel:

ClearAll[LL, CC, Dbuy, Dwait, Dlock, xx, xmin, xmax, n];
CC = p0^2/(p0 - LL);
Dbuy = (1 - CC)*(LL/p0) + (CC - p0)*Log[CC/p0];
Dwait = (1 - p1)*(LL/(p1 - p0));
Dlock = (1 - p1)*(1 - LL/p0 - LL/(p1 - p0)) + (p1 - CC)*(1 - LL/p0) - LL*Log[(p1 - p0)/(CC - p0)];

max[xx_, LLInit_, p0Init_, p1Init_] := 
  FindMaximum[{xx*(Dbuy*p0 + Dwait*0.4*p1 + Dlock*(p0 + 0.4*LL)),
    xx*(Dbuy + Dwait + Dlock) <= 1 && 0 <= LL <= p0 && p0 <= p1 <= 1}, 
  {{LL, LLInit}, {p0, p0Init}, {p1, p1Init}}]

xmin = 1.;
xmax = 6;
n = 20;
RevLocking1 = ConstantArray[{0, 0, 0, 0, 0}, n + 1];
Do[xx = xmin + (xmax - xmin) i/n;
 sol = If[i == 0, max[xx, 0.28, 0.95, 0.99],
   max[xx, RevLocking1[[i, 3]], RevLocking1[[i, 4]], RevLocking1[[i, 5]]]];
 RevLocking1[[i + 1]] = {xx, sol[[1]], LL, p0, p1} /. sol[[2]],
 {i, 0, n}]
RevLocking1 // TableForm

$$\left( \begin{array}{ccccc} 1. & 1.12066 & 0.373178 & 0.981281 & 1. \\ 1.25 & 1.10575 & 0.334785 & 0.979521 & 1. \\ 1.5 & 1.09425 & 0.304713 & 0.978577 & 1. \\ 1.75 & 1.0851 & 0.280353 & 0.978128 & 1. \\ 2. & 1.07763 & 0.260118 & 0.977986 & 1. \\ 2.25 & 1.07141 & 0.242977 & 0.97804 & 1. \\ 2.5 & 1.06615 & 0.228229 & 0.978219 & 1. \\ 2.75 & 1.06164 & 0.215377 & 0.978478 & 1. \\ 3. & 1.05773 & 0.204055 & 0.978788 & 1. \\ 3.25 & 1.0543 & 0.193992 & 0.979129 & 1. \\ 3.5 & 1.05127 & 0.184976 & 0.979488 & 1. \\ 3.75 & 1.04857 & 0.176844 & 0.979855 & 1. \\ 4. & 1.04614 & 0.169466 & 0.980226 & 1. \\ 4.25 & 1.04396 & 0.162735 & 0.980594 & 1. \\ 4.5 & 1.04198 & 0.156567 & 0.980958 & 1. \\ 4.75 & 1.04017 & 0.15089 & 0.981316 & 1. \\ 5. & 1.03852 & 0.145644 & 0.981665 & 1. \\ 5.25 & 1.037 & 0.140781 & 0.982006 & 1. \\ 5.5 & 1.0356 & 0.136258 & 0.982337 & 1. \\ 5.75 & 1.03431 & 0.132039 & 0.982659 & 1. \\ 6. & 1.03311 & 0.128093 & 0.982972 & 1. \\ \end{array} \right)$$

$\endgroup$
12
  • $\begingroup$ This is a promising start. I made a few changes in the maximization (there were some issues in the model) and it works just fine. Can we use NMaximize instead of FindMaximum though? $\endgroup$ – Ovunc Oct 20 '18 at 5:44
  • $\begingroup$ How are the models @TugrulTemel and I used different from your changed model? $\endgroup$ – JimB Oct 20 '18 at 13:25
  • $\begingroup$ I made a change in a constraint: Dbuy = (1 - CC)*(LL/p0) + (CC - p0) - p0*Log[CC/p0]; I also added a new constraint into the FindMaximum: p1*(p0 - LL) >= p0^2 Things work very good right now, and numbers make much more sense. I will go use my original method with these new constraints and let you know about it. Thanks so much! $\endgroup$ – Ovunc Oct 20 '18 at 16:53
  • $\begingroup$ Another change I made is "p0 + 0.4*LL" became "0.4*p0+LL", which is quite important:P $\endgroup$ – Ovunc Oct 20 '18 at 17:17
  • $\begingroup$ can you tell me how the following works? Do[xx = xmin + (xmax - xmin) i/n; sol = If[i == 0, max[xx, 0.28, 0.95, 0.99], max[xx, RevLocking1[[i, 3]], RevLocking1[[i, 4]], RevLocking1[[i, 5]]]]; RevLocking1[[i + 1]] = {xx, sol[[1]], LL, p0, p1} /. sol[[2]], {i, 0, n}] $\endgroup$ – Ovunc Oct 20 '18 at 17:53
1
$\begingroup$

You forgot to add x to the variable list. Also avoid using single capital letters as variable names since they may interfere with Mathematica built-in functions. I made some changes in your variable names, such x replaced with xx.

Try this:

ClearAll[LL, CC, Dbuy, Dwait, Dlock];
CC = p0^2/(p0 - LL) 
Dbuy = (1 - CC)*(LL/p0) + (CC - p0)*Log[CC/p0] 
Dwait = (1 - p1)*(LL/(p1 - p0)) 
Dlock = (1 - p1)*(1 - LL/p0 - LL/(p1 - p0)) + (p1 - CC)*(1 - LL/p0) - 
  LL*Log[(p1 - p0)/(CC - p0)]
RevLocking1 = 
 NMaximize[{xx*(Dbuy*p0 + Dwait*0.4*p1 + Dlock*(p0 + 0.4*LL)), 
   xx*(Dbuy + Dwait + Dlock) <= 1 && 0 <= LL <= p0 && 
    p0 <= p1 <= 1}, {p0, p1, LL, xx}]

giving the following result:

{0.147246, {LL -> 0.3322, p0 -> 0.571466, p1 -> 0.960381, 
  xx -> 0.457388}}

and it gives the following warning:

NMaximize::nrnum: "The function value -10.5461+1.01463\ I is not a real number at {LL,p0,p1,xx} = {0.21037,0.534614,0.529012,0.991907}. \!\(\*ButtonBox[\">>\",
Appearance->{Automatic, None},
BaseStyle->\"Link\",
ButtonData:>\"paclet:ref/NMaximize\",
ButtonNote->\"NMaximize::nrnum\"]\)"
$\endgroup$
2
  • $\begingroup$ x is not a variable, but a parameter. So I basically solve the same optimization problem for different values of x (between 1 and 6, with increments of 0.25). Thanks for the suggestions though, they are certainly helpful! $\endgroup$ – Ovunc Oct 20 '18 at 4:31
  • $\begingroup$ @Ovunc: You may want to try Manipulate for visualizing the solution with varying x. $\endgroup$ – Tugrul Temel Oct 20 '18 at 10:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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