# Avoid "Initial points not satisfying constraint"

I have this maximization problem which works for many examples I run before I made the For loops.

def[s_, cs_] := cs s
main[S_, cm_] := cm S
y[S_, Tp_, tech_] := Sqrt[S (Tp + tech)]

delta = 1; M = 2; V = 2; dec = 1; db = 10^3; ts = 0;
fs = 0; h = 0; re = 40^2; l = 0; p = 1; pm = 1; S = 1;
cb = 1; eps = 1; epsm = 1; tech = 1;

For[K = 1, K < 100, K += 50,
For[cs = 1, cs < 100, cs += 50,
For[cm = 1, cm < 100, cm += 50,
{tc1, s1, C11, b1} =
Values[Last[
NMaximize[
C1 - b db + h - cm (s + S) - re tc - cb (M + V) + epsm pm Sqrt[(s + dec S) (M + tech+V)],
Reduce[
{0 <= tc, tc <= M, C1 >= 0, s >= 0, b >= 0, K + eps p y[S, V + delta tc, tech] - def[s, cs] - main[S, cm] - C1 + b >= 0, C1 >= cb (M + V)},
{tc, s, C1, b}, Reals
],
{tc, s, C1, b}
]
]](*Values, Last*)
](*For*)
](*For*)
](*For*)


From my understanding satisfying

K + eps p y[S, V + delta tc, tech] - def[s, cs] - main[S, cm] - C1 + b >= 0


should not be a problem, because I allow b to be very high. So I really do not understand what the error message about "Initial Interval" is all about.

How do I have to change my code?

• Often a non-issue but can be: ? K reveals that K is a symbol used internally by Mathematica. Commented Feb 24, 2021 at 14:54
• At a minimum something is wrong with your syntax in NMaximize. Did you mean to include the equation and the Reduce expression in { }? Otherwise it just won't work at all. Furthermore, For returns Null, so you are continuously writing over the previous iteration's result and will only get the very last value. It has already been recommended to you to include Reap and Sow. Generally speaking, avoid loops unless you really really have to; you might use Table instead. Commented Feb 24, 2021 at 15:01

Block[{k = 51, cs = 51, cm = 51},
NMaximize[{C1 - b db + h - cm (s + S) - re tc - cb (M + V) +
epsm pm Sqrt[(s + dec S) (M + tech + V)],
Simplify[
And @@ {0 <= tc, tc <= M, C1 >= 0, s >= 0, b >= 0,
k + eps p y[S, V + delta tc, tech] - def[s, cs] - main[S, cm] -
C1 + b >= 0, C1 >= cb (M + V),
(s + dec S) (M + tech + V) >= 0, y[S, V + delta tc, tech]^2 >= 0}
]},
{tc, s, C1, b},
Method -> {"DifferentialEvolution", "InitialPoints" -> (
{tc, s, C1, b} /. FindInstance[{
k + eps p y[S, V + delta tc, tech] - def[s, cs] -
main[S, cm] - C1 + b >= 0,
(s + dec S) (M + tech + V) >= 0,
y[S, V + delta tc, tech]^2 >= 0},
{tc, s, C1, b},
20
])
},
MaxIterations -> 500
]
]

(*  {-2316.71, {tc -> 0., s -> 0., C1 -> 4., b -> 2.26795}}  *)


I had to add some missing {} and some constraints to keep the Sqrt terms real. The principal fix is to use FindInstance to find initial points.

• If you have an idea about bounds on the variables, you can use RandomPoint + ImplicitRegion, which seems a little faster: Block[{k = 51, cs = 51, cm = 51}, RandomPoint[ ImplicitRegion[ And @@ {k + eps p y[S, V + delta tc, tech] - def[s, cs] - main[S, cm] - C1 + b >= 0, (s + dec S) (M + tech + V) >= 0, y[S, V + delta tc, tech]^2 >= 0}, {tc, s, C1, b}], 20, {{0, 10}, {0, 10}, {0, 10}, {0, 10}}] ] Commented Feb 24, 2021 at 15:18
• Thank you very much. This worked out perfectly Commented Feb 26, 2021 at 8:57