# "Forcing output of a function to be Real in Mathematica" Non-linear optimization

I am trying to solve a Nonlinear optimization problem with Mathematica. Since the maximization problem includes exponentiation with fractional power, the result generates the complex numbers that stop the maximization. I used "Reduce" and "Rationalize" to get real values, but still, I get the following error.

NMaximize::nrnum: The function value "1.87153*10^9-105366. I" is not a real number at {y1,y2,y3,z1,z2,z3} = {0.102046,0.0181095,0.879844,0.273829,0.00536371,0.720807}.

Here is my code:

g=((500 y1+1595 y2-405 y3)^(0.4)*(700 z1+1757 z2-243 z3)^0.6*(0.4 (900 y1-121 y2+1879 y3)+0.6(900 z1-121 z2+1879 z3)))
Rationalize[g,0]
Reduce[(500 y1+1595 y2-405 y3)^(2/5) (700 z1+1757 z2-243 z3)^(3/5)==a && a∈Reals, Reals]

cons={500 y1+1595 y2-405 y3>=500 z1+1595 z2-405 z3  ,  700 z1+1757 z2-243 z3>=700 y1+1757 y2-243 y3    ,    y1+y2+y3==1 , z1+z2+z3==1 , y1>=0 , y2>=0 , y3>=0 , z1>=0 , z2>=0 , z3>=0}
cons3=Join[cons,{z1>=1/700 (-1757 z2+243 z3) && y1>=1/100 (-319 y2+81 y3)}]
vars={y1,y2,y3,z1,z2,z3}
sol3=Maximize[{g,cons3},vars]


I would appreciate if someone can help me.

• Why not just apply Maximize to the Re (real) part of your function? Aug 31, 2018 at 18:46
• ... and add Im[g]==0 to the constraints. Aug 31, 2018 at 19:15
• @David G. Stork thanks. that works now. I am just wondering about deleting Im[g]. Is it possible considering just Re[g], can change the optimal value? Aug 31, 2018 at 20:43
• Might want to redefine objective to be something very small if ordinary evaluation results in a non-real. Aug 31, 2018 at 21:38
• @Daniel Lichtblau could you please let me know how I can redefine objective?
– mah
Aug 31, 2018 at 21:51

Mathematica operates in the complex plane primarily, and thus fractional exponents typically result in complex answers. However, the Surd function allows you to specify that you are only interested in the real root of a number. Redefining g in terms of Surd rather than fractional powers immediately resolves the issue.

g = (Surd[500 y1 + 1595 y2 - 405 y3, 5]^2
* Surd[700 z1 + 1757 z2 - 243 z3, 5]^3
* (0.4 (900 y1 - 121 y2 + 1879 y3)
+ 0.6 (900 z1 - 121 z2 + 1879 z3)));

cons = {500 y1 + 1595 y2 - 405 y3 >= 500 z1 + 1595 z2 - 405 z3,
700 z1 + 1757 z2 - 243 z3 >= 700 y1 + 1757 y2 - 243 y3,
y1 + y2 + y3 == 1, z1 + z2 + z3 == 1, y1 >= 0, y2 >= 0, y3 >= 0,
z1 >= 0, z2 >= 0, z3 >= 0};
cons3 = Join[cons, {z1 >= 1/700 (-1757 z2 + 243 z3) && y1 >= 1/100 (-319 y2 + 81 y3)}];
vars = {y1, y2, y3, z1, z2, z3};
sol3 = Maximize[{g, cons3}, vars]


{613937., {y1 -> 1.807*10^-10, y2 -> 0.548779, y3 -> 0.451221, z1 -> 1.062*10^-9, z2 -> 0.548779, z3 -> 0.451221}}

Plugging this into your original g shows that this solution achieves the same numerical value of 613937, and is likely optimal.

gOrig = ((500 y1 + 1595 y2 - 405 y3)^(0.4)*(700 z1 + 1757 z2 - 243 z3)^0.6*(0.4 (900 y1 - 121 y2 + 1879 y3) + 0.6 (900 z1 - 121 z2 + 1879 z3)));
gOrig /. sol3[[2]]


613937.

I think that I've come across a similar problem before. I believe that Mathematica may need you to specify the variables in the domain of the function, otherwise it doesn't know the variables in respect of which it needs to perform optimisation.

When I modified your function as follows:

    g[y1_, y2_, y3_, z1_, z2_,
z3_] := ((500 y1 + 1595 y2 - 405 y3)^(0.4)*(700 z1 + 1757 z2 -
243 z3)^0.6*(0.4 (900 y1 - 121 y2 + 1879 y3) +
0.6 (900 z1 - 121 z2 + 1879 z3)))


and ran the rest of your code, I obtained the following result without error:

{g, {y1 -> 613048/853415, y2 -> 4617/1886000,
y3 -> 95311123/341366000, z1 -> 1160386913/1609540690,
z2 -> 2220777/1778498000, z3 -> 89428794763/321908138000}}


When I applied those values to g I obtained 410513.

By the way, I needed to restart the session after using your code. I did this using Quit[]