NMaximize and Conditional Loop

Question

I have developed the following code that simulates a process that decays over time and then returns to its initial state periodically.

When I pass a list of random integers in the interval $(1,0)$ as arguments, the function behaves as I expect it to.

However, when NMaximize calls the function, I get the following error:

NMaximize::nnum: The function value -1000 - 19 If[a[21] == 1,500, 1000] is not a number at {a[1], a[2], a[3], a[4], a[5], a[6], a[7], a[8], a[9], a[10], a[11], a[12], a[13], a[14], a[15], a[16], a[17], a[18], a[19], a[20]} = {1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}. >> The interesting part of the error is the complaint abouta[21]; this variable doesn't exist unless the function is called byNMaximize. The idea is thatNMaximize` will find integer, binary values for the parameters that will maximize the return value of the function.

My guess is that there is problem in how my function is being translated internally by NMaximize, but I am not not sure what to do at this point.

(*Assign Initial Values*)
Clear[f, i, a, vars, realconstraints, integerconstraints]
PeriodCapacityLoss = 10;
InitialCapacity = 1000;
OOSCapacity = 500;
AssymtoticCapacity = 200;
Periods = 20;
CurrentCapacity[1] = InitialCapacity;

(*Generate random series of cleaning flags*)
For[j = 1, j < Periods + 1, j++,
  RecoveryFlag[j] = RandomInteger[{0, 1}];
];

(*Function to simulate effect of capacity degredation and recovery*)
f[a_] := Module[{i},
  (* Set initial condition as clean *)  
  CurrentCapacity[1] = InitialCapacity;
  For[i = 2, i < Periods + 1, i++,
    CurrentCapacity[i] = 
      If[a[i] == 0 && a[i - 1] == 0, 
        CurrentCapacity[i - 1] - PeriodCapacityLoss, 
        If[a[i] == 1, OOSCapacity, InitialCapacity]
      ];
  ];
  Return[Total[Map[CurrentCapacity, Range[Periods]]]];
];

(*Pass random cleaning flags to degredation function and plot*)
Print[f[RecoveryFlag]];

ListLinePlot[Array[CurrentCapacity, Periods]]

(* Define the integer constraints *)
vars = Array[a, Periods];
realconstraints = And @@ Map[Greater[2, #, 0]&, vars];
integerconstraints = Append[realconstraints, Element[vars, Integers]];

(* 
  Find the value of the recovery flag that maximizes capacity across the 
  time window
*)
NMaximize[{f[a], integerconstraints}, vars, Method->{"DifferentialEvolution"}]

Answer 1

Update

Aha, 8 years passed and now I'm able to figure out what's wrong with OP's original attempt. The problem can be boiled down to the following:

Why does an a[21] appear in the warning message?

This is actually a matter of evaluation order. Just observe the output of the following sample:

Clear[a, i];
i = 1;
If[a[i] > 0, 1, If[a[i] < 2, 4]]
(* If[a[1] > 0, 1, If[a[i] < 2, 4]] *)

As we can see, there's an a[i] in 2nd If[…], this is because If owns HoldRest attribute, so the a[i] inside 2nd If[…] is never evaluated.

"So what? This won't cause any problem, because a[i] will finally evaluate to a[1] at sometime! " Yes, if only the value of i never changes afterwards, but sadly it's not the case for OP's f: the value of i is changing in the For loop, which finally becomes 21!

So, how to fix? Using PiecewiseExpand as shown below is of course a solution, but a more on-target solution is to adjust the evaluation order to make the a[i]s evaluate at proper timings, which can be done with With:

f[a_] := Module[{i},(*Set initial condition as clean*)
  CurrentCapacity[1] = InitialCapacity;
  For[i = 2, i < Periods + 1, i++, 
   CurrentCapacity[i] =(*PiecewiseExpand@*)
     With[{i = i}, 
      If[a[i] == 0 && a[i - 1] == 0, CurrentCapacity[i - 1] - PeriodCapacityLoss, 
       If[a[i] == 1, OOSCapacity, InitialCapacity]]];];
  Return[Total[Map[CurrentCapacity, Range[Periods]]]];]

You'll see some of the i becomes red, it's merely for warning. If you don't like it, just change the With[…] to With[{ii = i}, If[a[ii] == 0 && a[ii - 1] == 0, CurrentCapacity[ii - 1] - PeriodCapacityLoss, If[a[ii] == 1, OOSCapacity, InitialCapacity]]].

---

Maybe there're deeper reasons, but I can't figure it out right now: it seems that NMinimize doesn't like If. Add a PiecewiseExpand to the definition of f will fix the problem:

f[a_] := Module[{i},(*Set initial condition as clean*)
   CurrentCapacity[1] = InitialCapacity;
   For[i = 2, i < Periods + 1, i++, 
    CurrentCapacity[i] = 
      PiecewiseExpand@
       If[a[i] == 0 && a[i - 1] == 0, 
        CurrentCapacity[i - 1] - PeriodCapacityLoss, 
        If[a[i] == 1, OOSCapacity, InitialCapacity]];];
   Return[Total[Map[CurrentCapacity, Range[Periods]]]];];

vars = Array[a, Periods];
realconstraints = And @@ Map[Greater[2, #, 0] &, vars];
integerconstraints = Append[realconstraints, Element[vars, Integers]];

NMaximize[{f[a], integerconstraints}, vars, Method -> {"DifferentialEvolution"}]

{10500., {a[1] -> 1, a[2] -> 1, a[3] -> 1, a[4] -> 1, a[5] -> 1,
a[6] -> 1, a[7] -> 1, a[8] -> 1, a[9] -> 1, a[10] -> 1, a[11] -> 1,
a[12] -> 1, a[13] -> 1, a[14] -> 1, a[15] -> 1, a[16] -> 1, a[17] -> 1, a[18] -> 1, a[19] -> 1, a[20] -> 1}}