# “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? – David G. Stork Aug 31 '18 at 18:46
• ... and add Im[g]==0 to the constraints. – Henrik Schumacher Aug 31 '18 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? – Soodi Aug 31 '18 at 20:43
• Might want to redefine objective to be something very small if ordinary evaluation results in a non-real. – Daniel Lichtblau Aug 31 '18 at 21:38
• @Daniel Lichtblau could you please let me know how I can redefine objective? – mah Aug 31 '18 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[]


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