Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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"}]
share|improve this question
    
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] := ... –  Daniel Lichtblau Aug 1 '13 at 5:50
    
@DanielLichtblau But in this sample, a is only a Head –  xzczd Aug 1 '13 at 5:54
    
@xzczd I had missed that. Would be better to make it a vector, say, or an explicit function. –  Daniel Lichtblau Aug 1 '13 at 18:04

1 Answer 1

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

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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