Coverage for moptipyapps / prodsched / rop_experiment.py: 33%

1""" 

2A small template for ROP-based experiments. 

3 

4This experiment uses the NSGA-II algorithm to optimize Re-Order-Points (ROPs) 

5to achieve both a high worst-case fillrate and a low worst-case average stock 

6level. 

7 

8It is just a preliminary idea. 

9 

10In manufacturing systems, the concept of re-order points (ROPs) can be used 

11to schedule production. 

12The ROP basically provides one limit `X` per product. 

13If there are less or equal (`<=`) `X` units of the product in the warehouse, 

14a production order is issued so that one new unit is produced. 

15 

16The question is how should `X` be set so that we can 

17 

181. satisfy as many customer demands as possible directly when they come in, 

19 i.e., maximize the fill rate 

20 (:attr:`~moptipyapps.prodsched.statistics.Statistics.immediate_rates`, 

21 represented by 

22 :mod:`~moptipyapps.prodsched.objectives.worst_and_mean_fill_rate`) and 

232. have a low total average number of product units sitting around in our 

24 warehouse, i.e., minimize the stock level 

25 (:attr:`~moptipyapps.prodsched.statistics.Statistics.stock_levels`, 

26 represented by 

27 :mod:`~moptipyapps.prodsched.objectives.max_stocklevel`). 

28 

29Since we have `n=10` products in the Thürer [1] base scenario, we also have 

30`10` such ROP values `X`. 

31 

32Now how can we know the fill rate and the stock level? 

33This is done by simulating the whole system in action. 

34We therefore use instances generated using the Thürer-style distributions as 

35defined in :mod:`~moptipyapps.prodsched.mfc_generator` [1]. 

36 

37Our simulations (see :mod:`~moptipyapps.prodsched.simulation`) are not 

38randomized. 

39Instead, each of them is based on a fixed 

40:mod:`~moptipyapps.prodsched.instance`. 

41Each :mod:`~moptipyapps.prodsched.instance.Instance` defines exactly when a 

42customer :class:`~moptipyapps.prodsched.instance.Demand` comes in and, 

43for many time windows, the production time needed by a machine to produce 

44one unit of product. 

45Of course, all of these values follow the random distributions (Erlang, Gamma, 

46with respective parameters) given in Tables 2 and 3 of the original Thürer 

47paper [1] and implemented in :mod:`~moptipyapps.prodsched.mfc_generator`. 

48But apart from this, they are fixed per instance. 

49 

50This means that we can take the same ROP and run the simulation twice for a 

51given instance and will get exactly the same results and 

52:mod:`~moptipyapps.prodsched.statistics` (fill rates, stock levels, etc.). 

53 

54Of course, using a single fixed :mod:`~moptipyapps.prodsched.instance` may 

55be misleading. 

56Maybe we would think that a certain ROP is very good ... but it is only 

57good on the specific instance we tried. 

58 

59So as a second step, we here generate 11 instances (via 

60:func:`~moptipyapps.prodsched.instances.get_instances`). 

61And then we look at the worst fill rate and the worst stock level over all 11 

62instances (more or less). 

63And we use that to judge whether an ROP is good or not. 

64 

65This leaves the question: 

66Where do these ROPs come from? 

67 

68They come from an optimization process. 

69First, we define that ROPs be integer vectors (i.e., from an 

70:class:`~moptipy.spaces.intspace.IntSpace`) where each element comes from 

71the range `0..63`. 

72We sample the initial solutions randomly from that interval 

73(via :class:`~moptipy.operators.intspace.op0_random.Op0Random`). 

74 

75As unary search operator, we take an existing ROP and, for a number of 

76elements, sample a new value normally distributed around it (with standard 

77deviation 2.5) but rounded to integer. 

78We do this for a binomially distributed number of elements, exactly like the 

79bit-string based (1+1) EA would do it, but implemented for integers by 

80operator :class:`~moptipy.operators.intspace.op1_mnormal.Op1MNormal`. 

81 

82As binary operator, we use uniform crossover, given as 

83:class:`~moptipy.operators.intspace.op2_uniform.Op2Uniform`. 

84This operator takes two existing solutions and creates a new solution by 

85element-wise copying elements of the parents. 

86For each element, it randomly decides from which parent solution it should be 

87copied. 

88 

89As optimization algorithm, we use NSGA-II implemented by class 

90:class:`~moptipy.algorithms.mo.nsga2.NSGA2`. 

91This is a multi-objective optimization algorithm. 

92We do this because we have two goals: 

93 

941. Maximize the worst-case fill rates, 

952. Minimize the worst-case average stock level. 

96 

97Regarding the first objective, we have a small tweak: 

98Assume that, over all 11 instances, `PM` be the worst fill rate per product 

99(in [0,1], 0 being worst) and `AM` be the worst average fill rate over all 

100products (0 worst, 1 best). 

101Then our objective value -- subject to minimization -- is 

102`(1 - PM) * 100 + (1 - AM)`. 

103This is implemented in module 

104:mod:`~moptipyapps.prodsched.objectives.worst_and_mean_fill_rate`. 

105 

106The objective function minimizing the stock level is implemented in module 

107:mod:`~moptipyapps.prodsched.objectives.max_stocklevel`. 

108 

109In summary, what we do is this: 

110 

1111. The optimization algorithm proposes ROPs, each of which being an integer 

112 vector with 10 values (1 value per product). 

113 (:class:`~moptipy.spaces.intspace.IntSpace`) 

114 These integer vectors are the elements of the search space. 

115 

1162. The ROP is evaluated by simulating it 11 times (using 10000 time units per 

117 instance, 3000 of which are used for warmup). 

118 This is implemented as an :mod:`~moptipy.api.encoding`, 

119 :class:`~moptipyapps.prodsched.rop_multisimulation.ROPMultiSimulation`. 

120 

1213. This :mod:`~moptipyapps.prodsched.simulation`-based decoding procedure maps 

122 the ROP vectors to the solution space. In our case, this solution space are 

123 just multi-statistics records, as implemented in module 

124 :mod:`~moptipyapps.prodsched.multistatistics`, where a corresponding 

125 :mod:`~moptipy.api.space` implementation is also provided. 

126 

1274. Each of the 11 simulations is based on one fixed 

128 :mod:`~moptipyapps.prodsched.instance`, where all customer demands and 

129 machine work times (at certain time intervals) are fixed (but were sampled 

130 based on the distributions given in the Thürer paper [1]) using module 

131 :mod:`~moptipyapps.prodsched.mfc_generator`. 

132 

1335. Each of the 11 simulations has per-product fill rates `PM_i,j`, one average 

134 fill rate `AM_i`, one average stock level `SL_i` stored in the 

135 :mod:`~moptipyapps.prodsched.multistatistics` records. 

136 

1376. As first objective, we use the smallest `PM_i_j` as `PM` and the smallest 

138 `AM_i` as `AM` and compute `(1 - AM) * 100 + (1 - PM)` in 

139 :mod:`~moptipyapps.prodsched.objectives.worst_and_mean_fill_rate`. 

140 

1417. As second objective, we use the largest `SL_i` in 

142 :mod:`~moptipyapps.prodsched.objectives.max_stocklevel`. 

143 

1448. NSGA-II, given in :mod:`~moptipy.algorithms.mo.nsga2`, maintains a 

145 population of solutions which it evaluates like that. 

146 

1479. NSGA-II decides which solutions to keep based on the current Pareto front 

148 and crowding distance. 

149 

15010. The retained solutions are reproduced, either via unary or binary search 

151 operators. 

152 

15311. The unary search operator changes a random number (binomial distribution) 

154 of elements of a ROP vector (via normal distribution). 

155 It is given in 

156 :class:`~moptipy.operators.intspace.op1_mnormal.Op1MNormal`. 

157 

15812. The binary search operator is simple uniform crossover, i.e., fills a 

159 new solution with elements of either of the two parent solutions. It is 

160 defined in :class:`~moptipy.operators.intspace.op2_uniform.Op2Uniform`. 

161 

162The core idea is that we do not use randomized ARENA-like simulation. 

163Instead, we use a simulator that is based on deterministic, fully pre-defined 

164instances. 

165These instances are still randomly generated according to the distributions 

166given by Thürer in [1]. 

167However, they have all values pre-determined. 

168This allows us to run a simulation with the same ROP twice and get the exactly 

169same result. 

170If two different ROPs are simulated, but both of them decide to produce 1 unit 

171of product `A` on machine `V` at time unit `T`, then for both of them this 

172will take exactly the same amount of time. 

173 

174Simulation results are thus less noisy. 

175 

176Of course, there is a danger of overfitting. 

177This is why we need to use multiple simulations to check a ROP. 

178And then we take the worst-case results. 

179 

180So this is the idea. 

181 

1821. Matthias Thürer, Nuno O. Fernandes, Hermann Lödding, and Mark Stevenson. 

183 Material Flow Control in Make-to-Stock Production Systems: An Assessment of 

184 Order Generation, Order Release and Production Authorization by Simulation 

185 Flexible Services and Manufacturing Journal. 37(1):1-37. March 2025. 

186 doi: https://doi.org/10.1007/s10696-024-09532-2 

187""" 

188 

189 

190import argparse 

191from typing import Final 

192 

193from moptipy.algorithms.mo.nsga2 import NSGA2 

194from moptipy.api.experiment import run_experiment 

195from moptipy.api.mo_execution import MOExecution 

196from moptipy.mo.problem.weighted_sum import WeightedSum 

197from moptipy.operators.intspace.op0_random import Op0Random 

198from moptipy.operators.intspace.op1_mnormal import Op1MNormal 

199from moptipy.operators.intspace.op2_uniform import Op2Uniform 

200from moptipy.spaces.intspace import IntSpace 

201from pycommons.io.console import logger 

202from pycommons.io.path import Path 

203from pycommons.types import check_int_range 

204 

205from moptipyapps.prodsched.instance import Instance 

206from moptipyapps.prodsched.instances import get_instances 

207from moptipyapps.prodsched.multistatistics import MultiStatisticsSpace 

208from moptipyapps.prodsched.objectives.max_stocklevel import MaxStockLevel 

209from moptipyapps.prodsched.objectives.worst_and_mean_fill_rate import ( 

210 WorstAndMeanFillRate, 

211) 

212from moptipyapps.prodsched.rop_multisimulation import ROPMultiSimulation 

213from moptipyapps.utils.shared import moptipyapps_argparser 

214 

215 

216def run(dest: str, instances: str, n_inst: int, n_runs: int, 

217 max_fes: int, ps: int) -> None: 

218 """ 

219 Run the experiment. 

220 

221 :param dest: the destination directory 

222 :param instances: the directory with the instances 

223 :param n_inst: the number of instances 

224 :param n_runs: the number of runs 

225 :param max_fes: the maximum FEs 

226 :param ps: the population size 

227 """ 

228 logger(f"Beginning experiment with dest={dest!r}, instances={instances!r}" 

229 f", n_inst={n_inst}, n_runs={n_runs}, and max_fes={max_fes}.") 

230 use_dest: Final[Path] = Path(dest) 

231 use_dest.ensure_dir_exists() 

232 logger(f"Destination folder is {use_dest!r}.") 

233 

234 use_insts: Final[Path] = Path(instances) 

235 use_insts.ensure_dir_exists() 

236 logger(f"Instances folder is {use_insts!r}.") 

237 

238 check_int_range(n_inst, "n_inst", 1, 128) 

239 check_int_range(max_fes, "max_fes", 10, 10 ** 10) 

240 check_int_range(ps, "ps", 4, 16384) 

241 

242 logger(f"Loading {n_inst} instances from {use_insts!r}.") 

243 insts: Final[tuple[Instance, ...]] = get_instances(n_inst, instances) 

244 if tuple.__len__(insts) != n_inst: 

245 raise ValueError("Could not load required instances.") 

246 logger(f"Loaded {n_inst} instances from {use_insts!r}.") 

247 

248 n_prod: int | None = None 

249 for inst in insts: 

250 if n_prod is None: 

251 n_prod = inst.n_products 

252 elif n_prod != inst.n_products: 

253 raise ValueError("Inconsistent number of products!") 

254 if n_prod is None: 

255 raise ValueError("No instances?") 

256 

257 search_space: Final[IntSpace] = IntSpace(n_prod, 0, 63) 

258 op0: Final[Op0Random] = Op0Random(search_space) 

259 op1: Final[Op1MNormal] = Op1MNormal(search_space, sd=2.5) 

260 op2: Final = Op2Uniform() 

261 algo: Final[NSGA2] = NSGA2(op0, op1, op2, ps, 1 / min(16, ps)) 

262 solution_space: Final[MultiStatisticsSpace] = MultiStatisticsSpace(insts) 

263 encoding: Final[ROPMultiSimulation] = ROPMultiSimulation(solution_space) 

264 f1: Final[WorstAndMeanFillRate] = WorstAndMeanFillRate() 

265 f2: Final[MaxStockLevel] = MaxStockLevel() 

266 ws: Final[WeightedSum] = WeightedSum((f1, f2), ( 

267 (1 / (f1.upper_bound() - f1.lower_bound())), 1 / (2 * n_prod))) 

268 

269 def __setup(_) -> MOExecution: 

270 """ 

271 Set up the experiment. 

272 

273 :return: the execution 

274 """ 

275 return (MOExecution() 

276 .set_search_space(search_space) 

277 .set_algorithm(algo) 

278 .set_solution_space(solution_space) 

279 .set_objective(ws) 

280 .set_encoding(encoding) 

281 .set_max_fes(max_fes) 

282 .set_log_improvements(True)) 

283 

284 run_experiment(base_dir=use_dest, instances=(lambda: "all", ), 

285 setups=(__setup, ), n_runs=n_runs, 

286 pre_warmup_fes=2, perform_warmup=False, 

287 perform_pre_warmup=True) 

288 

289 

290# Run the experiment from the command line 

291if __name__ == "__main__": 

292 parser: Final[argparse.ArgumentParser] = moptipyapps_argparser( 

293 __file__, "ROP-based MFC Optimization", 

294 "Run a small experiment with ROP-based MFC optimization.") 

295 parser.add_argument( 

296 "dest", help="the directory to store the experimental results under", 

297 type=Path, nargs="?", default="./results/") 

298 parser.add_argument( 

299 "insts", help="the directory with the instances", 

300 type=Path, nargs="?", default="./instances/") 

301 parser.add_argument( 

302 "n_inst", help="the number of instances", 

303 type=int, nargs="?", default=11) 

304 parser.add_argument( 

305 "n_runs", help="the number of runs", 

306 type=int, nargs="?", default=31) 

307 parser.add_argument( 

308 "max_fes", help="the number of FEs per run", 

309 type=int, nargs="?", default=8192) 

310 parser.add_argument( 

311 "ps", help="the population size of NSGA-II", 

312 type=int, nargs="?", default=64) 

313 args: Final[argparse.Namespace] = parser.parse_args() 

314 run(args.dest, args.insts, args.n_inst, args.n_runs, args.max_fes, 

315 args.ps)