2
$\begingroup$

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"}]
$\endgroup$
3
  • $\begingroup$ Possibly NMaximize will be less inclined to sulk if the function is a "black box", that is, only evaluating for explicitly numeric arguments. That can be attained via e.g. f[a_?NumberQ] := ... $\endgroup$ Aug 1, 2013 at 5:50
  • $\begingroup$ @DanielLichtblau But in this sample, a is only a Head $\endgroup$
    – xzczd
    Aug 1, 2013 at 5:54
  • $\begingroup$ @xzczd I had missed that. Would be better to make it a vector, say, or an explicit function. $\endgroup$ Aug 1, 2013 at 18:04

1 Answer 1

1
$\begingroup$

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}}

$\endgroup$
5
  • $\begingroup$ All PiecwiseExpand seems to be doing here is defining a default value so the If statement doesn't remain unevaluated; e.g. providing 0 as the last argument of the last If gives the same fix. $\endgroup$ Oct 7, 2021 at 10:45
  • $\begingroup$ @RonaldMonson Not the same. The maximum in my case is 10500, but yours is 1000. BTW I manage to figure out the root of OP's problem, check my edit. $\endgroup$
    – xzczd
    Oct 7, 2021 at 11:28
  • $\begingroup$ I'm glad this seems to be sorted - I haven't gone through the details (my eyes glaze over when imperative algorithms are used these days - and indeed such bugs are an illustration of the superiority of functional solutions IMO). But the broader point that led to this re-visiting remains: PiecewiseExpand's transformations weren't strictly necessary in the light of a more straightforward re-factoring. $\endgroup$ Oct 8, 2021 at 0:55
  • $\begingroup$ @RonaldMonson I'd say it's natural that PiecewiseExpand isn't necessary. (It doesn't even exist in programming languages like C, right? :D ) But at least in this case it provides a sweet alternative for fixing the problem. (PiecewiseExpand@ is easier to type than With[{i = i}, …] in my opinion. ) $\endgroup$
    – xzczd
    Oct 8, 2021 at 1:35
  • $\begingroup$ I'd say it is problematic when PiecewiseExpand is used unnecessarily. It is potentially much more expensive given the rich transformation rules it contains. Don't get me wrong, I think it can be a very powerful tool in mathematical settings (unlike other languages, indeed) and my original question was exploring if it could be used likewise beyond such settings. The issue seem to be that it is performing a dual role when applied to procedural expressions, casting them into mathematical form and performing simplifications/expansions. Sometimes it is only the former that is intended. $\endgroup$ Oct 8, 2021 at 2:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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