Coverage for agentlib_flexquant/modules/shadow_mpc.py: 71%

1""" 

2Defines shadow MPC and MINLP-MPC for positive/negative flexibility quantification. 

3""" 

4import os 

5import math 

6import numpy as np 

7import pandas as pd 

8from pydantic import Field 

9from typing import Dict, Union, Optional 

10from collections.abc import Iterable 

11from agentlib.core.datamodels import AgentVariable, Source 

12from agentlib_mpc.modules.mpc import mpc_full, minlp_mpc 

13from agentlib_flexquant.utils.data_handling import fill_nans, MEAN 

14from agentlib_flexquant.data_structures.globals import (full_trajectory_suffix, 

15 base_vars_to_communicate_suffix) 

16import agentlib_flexquant.data_structures.globals as glbs 

17from agentlib_flexquant.optimization_backends.constrained_cia import ConstrainedCasADiCIABackend 

18 

19 

20class FlexibilityShadowMPCConfig(mpc_full.MPCConfig): 

21 

22 baseline_input_names: list[str] = Field(default=[]) 

23 custom_input_names: list[Dict] = Field(default=[]) 

24 baseline_state_names: list[str] = Field(default=[]) 

25 full_control_names: list[str] = Field(default=[]) 

26 

27 

28 baseline_agent_id: str = "" 

29 

30 casadi_sim_time_step: int = Field( 

31 default=0, 

32 description="Time step for simulation with Casadi simulator. " 

33 "Value is read from FlexQuantConfig", 

34 ) 

35 power_variable_name: str = Field( 

36 default=None, description="Name of the power variable in the " 

37 "shadow mpc model." 

38 ) 

39 storage_variable_name: Optional[str] = Field( 

40 default=None, description="Name of the storage variable in the " 

41 "shadow mpc model." 

42 ) 

43 

44 

45class FlexibilityShadowMPC(mpc_full.MPC): 

46 """Shadow MPC for calculating positive/negative flexibility offers.""" 

47 

48 config: FlexibilityShadowMPCConfig 

49 

50 def __init__(self, *args, **kwargs): 

51 # initialize flex_results with None 

52 self.flex_results = None 

53 

54 super().__init__(*args, **kwargs) 

55 

56 # setup look up dict to track incoming inputs and states 

57 # (maps name as str to actual AgentVariable) 

58 input_names_list = [var["name"] for var in self.config.custom_input_names] 

59 self._track_base_comm_vars_dict: Dict[str, Union[AgentVariable, None]] = {} 

60 for comm_var in self.config.inputs + self.config.states: 

61 if (comm_var.name in self.config.full_control_names or 

62 comm_var.name + base_vars_to_communicate_suffix in 

63 self.config.baseline_input_names or 

64 comm_var.name + base_vars_to_communicate_suffix in 

65 self.config.baseline_state_names or 

66 comm_var.name in input_names_list): 

67 comm_var.value = None 

68 self._track_base_comm_vars_dict[comm_var.name] = comm_var.copy(deep=True) 

69 # set up necessary components if simulation is enabled 

70 if self.config.casadi_sim_time_step > 0: 

71 # generate a separate simulation model for integration to ensure 

72 # the model used in MPC optimization remains unaffected 

73 self.flex_model = type(self.model)(dt=self.config.casadi_sim_time_step) 

74 # generate the filename for the simulation results 

75 self.res_file_flex = self.config.optimization_backend["results_file"].replace( 

76 "_flex", "_sim_flex" 

77 ) 

78 # clear the casadi simulator result at the first time step if already exists 

79 try: 

80 os.remove(self.res_file_flex) 

81 except FileNotFoundError: 

82 pass 

83 

84 def set_output(self, solution): 

85 """Takes the solution from optimization backend and sends it to AgentVariables.""" 

86 # Output must be defined in the config as "type"="pd.Series" 

87 if not self.config.set_outputs: 

88 return 

89 self.logger.info("Sending optimal output values to data_broker.") 

90 df = solution.df 

91 self.sim_flex_model(solution) 

92 if self.flex_results is not None: 

93 for output in self.var_ref.outputs: 

94 if output not in [ 

95 self.config.power_variable_name, 

96 self.config.storage_variable_name, 

97 ]: 

98 series = df.variable[output] 

99 self.set(output, series) 

100 # send the power and storage variable value from simulation results 

101 upsampled_output_power = self.flex_results[self.config.power_variable_name] 

102 self.set(self.config.power_variable_name, upsampled_output_power) 

103 if self.config.storage_variable_name is not None: 

104 upsampled_output_storage = self.flex_results[self.config.storage_variable_name] 

105 self.set(self.config.storage_variable_name, upsampled_output_storage.dropna()) 

106 else: 

107 for output in self.var_ref.outputs: 

108 series = df.variable[output] 

109 self.set(output, series) 

110 

111 def sim_flex_model(self, solution): 

112 """simulate the flex model over the preditcion horizon and save results""" 

113 

114 # return if sim_time_step is not a positive integer and system is in provision 

115 if not (self.config.casadi_sim_time_step > 0 and not self.get(glbs.PROVISION_VAR_NAME).value): 

116 return 

117 

118 # read the defined simulation time step 

119 sim_time_step = self.config.casadi_sim_time_step 

120 mpc_time_step = self.config.time_step 

121 

122 # set the horizon length and the number of simulation steps 

123 total_horizon_time = int(self.config.prediction_horizon * self.config.time_step) 

124 n_simulation_steps = math.ceil(total_horizon_time / sim_time_step) 

125 

126 # read the current optimization result 

127 result_df = solution.df 

128 

129 # initialize the flex sim results Dataframe 

130 self._initialize_flex_results( 

131 n_simulation_steps, total_horizon_time, sim_time_step, result_df 

132 ) 

133 

134 # Update model parameters and initial states 

135 self._update_model_parameters() 

136 self._update_initial_states(result_df) 

137 

138 # Run simulation 

139 self._run_simulation( 

140 n_simulation_steps, sim_time_step, mpc_time_step, result_df, total_horizon_time 

141 ) 

142 

143 # set index of flex results to the same as mpc result 

144 store_results_df = self.flex_results.copy(deep=True) 

145 store_results_df.index = self.flex_results.index.tolist() 

146 

147 # save results 

148 if not os.path.exists(self.res_file_flex): 

149 store_results_df.to_csv(self.res_file_flex) 

150 else: 

151 store_results_df.to_csv(self.res_file_flex, mode="a", header=False) 

152 

153 # set the flex results format same as mpc result while updating Agentvariable 

154 self.flex_results.index = self.flex_results.index.get_level_values(1) 

155 

156 def register_callbacks(self): 

157 for control_var in self.config.controls: 

158 self.agent.data_broker.register_callback( 

159 name=control_var.name + full_trajectory_suffix, 

160 alias=control_var.name + full_trajectory_suffix, 

161 callback=self.calc_flex_callback, 

162 source=Source(agent_id=self.config.baseline_agent_id, module_id=None) 

163 ) 

164 for base_inputs in self.config.baseline_input_names: 

165 self.agent.data_broker.register_callback( 

166 name=base_inputs.removesuffix(base_vars_to_communicate_suffix), # update MPC variable 

167 alias=base_inputs, 

168 callback=self.calc_flex_callback, 

169 source=Source(agent_id=self.config.baseline_agent_id, module_id=None) 

170 ) 

171 for custom_inputs in self.config.custom_input_names: 

172 self.agent.data_broker.register_callback( 

173 name=custom_inputs["name"], 

174 alias=custom_inputs["alias"], 

175 callback=self.calc_flex_callback 

176 ) 

177 for base_states in self.config.baseline_state_names: 

178 self.agent.data_broker.register_callback( 

179 name=base_states.removesuffix(base_vars_to_communicate_suffix), # update MPC variable 

180 alias=base_states, 

181 callback=self.calc_flex_callback, 

182 source=Source(agent_id=self.config.baseline_agent_id, module_id=None) 

183 ) 

184 super().register_callbacks() 

185 

186 def calc_flex_callback(self, inp: AgentVariable, name: str): 

187 """Ensure that all control trajectories and Baseline inputs/states 

188 have been set before starting the calculation. 

189 

190 """ 

191 # during provision do not calculate flex 

192 if self.get(glbs.PROVISION_VAR_NAME).value: 

193 return 

194 

195 # do not trigger callback on self set variables 

196 if self.agent.config.id == inp.source.agent_id: 

197 return 

198 

199 # get the value of the input 

200 vals = inp.value 

201 

202 if inp.name in self.config.full_control_names: 

203 if vals.isna().any(): 

204 vals = fill_nans(series=vals, method=MEAN) 

205 # add time shift env.time to the incoming variable to adapt to mpc output, 

206 # which starts at t=0 

207 if vals.index[0] == 0: 

208 self.logger.info(f"The incoming variable {inp.name} starts with a time " 

209 f"index of 0. Adding the current environment time.") 

210 vals.index += self.env.time 

211 

212 # update value in the tracking dictionary 

213 self._track_base_comm_vars_dict[name].value = vals 

214 # set value 

215 self.set(name, vals) 

216 

217 # make sure all necessary inputs are set 

218 if all(x.value is not None for x in self._track_base_comm_vars_dict.values()): 

219 self.do_step() 

220 for _, comm_var in self._track_base_comm_vars_dict.items(): 

221 comm_var.value = None 

222 

223 def process(self): 

224 # the shadow mpc should only be run after the results of the baseline are sent 

225 yield self.env.event() 

226 

227 def _initialize_flex_results( 

228 self, n_simulation_steps, horizon_length, sim_time_step, result_df 

229 ): 

230 """Initialize the flex results dataframe with the correct dimension 

231 and index and fill with existing results from optimization 

232 

233 """ 

234 

235 # create MultiIndex for collocation points 

236 index_coll = pd.MultiIndex.from_arrays( 

237 [[self.env.now] * len(result_df.index), result_df.index], 

238 names=["time_step", "time"] 

239 # Match the names with multi_index but note they're reversed 

240 ) 

241 # create Multiindex for full simulation sample times 

242 index_full_sample = pd.MultiIndex.from_tuples( 

243 zip( 

244 [self.env.now] * (n_simulation_steps + 1), 

245 range(0, horizon_length + sim_time_step, sim_time_step), 

246 ), 

247 names=["time_step", "time"], 

248 ) 

249 # merge indexes 

250 new_index = index_coll.union(index_full_sample).sort_values() 

251 # initialize the flex results with correct dimension 

252 self.flex_results = pd.DataFrame(np.nan, 

253 index=new_index, 

254 columns=self.var_ref.outputs) 

255 

256 # Get the optimization outputs and create a series for fixed 

257 # optimization outputs with the correct MultiIndex format 

258 opti_outputs = result_df.variable[self.config.power_variable_name] 

259 fixed_opti_output = pd.Series( 

260 opti_outputs.values, 

261 index=index_coll, 

262 ) 

263 # fill the output value at the time step where it already exists 

264 # in optimization output 

265 for idx in fixed_opti_output.index: 

266 if idx in self.flex_results.index: 

267 self.flex_results.loc[idx, self.config.power_variable_name] = ( 

268 fixed_opti_output)[idx] 

269 

270 def _update_model_parameters(self): 

271 """update the value of module parameters with value from config, 

272 since creating a model just reads the value in the model class 

273 but not the config. 

274 

275 """ 

276 

277 for par in self.config.parameters: 

278 self.flex_model.set(par.name, par.value) 

279 

280 def _update_initial_states(self, result_df): 

281 """set the initial value of states""" 

282 

283 # get state values from the mpc optimization result 

284 state_values = result_df.variable[self.var_ref.states] 

285 # update state values with last measurement 

286 for state, value in zip(self.var_ref.states, state_values.iloc[0]): 

287 self.flex_model.set(state, value) 

288 

289 def _run_simulation( 

290 self, n_simulation_steps, sim_time_step, mpc_time_step, result_df, total_horizon_time 

291 ): 

292 """simulate with flex model over the prediction horizon 

293 

294 """ 

295 

296 # get control and input values from the mpc optimization result 

297 control_values = result_df.variable[self.var_ref.controls].dropna() 

298 input_values = result_df.parameter[self.var_ref.inputs].dropna() 

299 

300 # Get the simulation time step index 

301 sim_time_index = np.arange(0, (n_simulation_steps + 1) * sim_time_step, sim_time_step) 

302 

303 # Reindex the controls and inputs to sim_time_index 

304 control_values_full = control_values.copy().reindex(sim_time_index, method="ffill") 

305 input_values_full = input_values.copy().reindex(sim_time_index, method="nearest") 

306 

307 for i in range(0, n_simulation_steps): 

308 current_sim_time = i * sim_time_step 

309 

310 # Apply control and input values from the appropriate MPC step 

311 for control, value in zip( 

312 self.var_ref.controls, control_values_full.loc[current_sim_time] 

313 ): 

314 self.flex_model.set(control, value) 

315 

316 for input_var, value in zip( 

317 self.var_ref.inputs, input_values_full.loc[current_sim_time] 

318 ): 

319 # change the type of iterable input, since casadi model can't deal with iterable 

320 if issubclass(eval(self.flex_model.get(input_var).type), Iterable): 

321 self.flex_model.get(input_var).type = type(value).__name__ 

322 self.flex_model.set(input_var, value) 

323 

324 # do integration 

325 # reduce the simulation time step so that the total horizon time will not be exceeded 

326 if current_sim_time + sim_time_step <= total_horizon_time: 

327 t_sample = sim_time_step 

328 else: 

329 t_sample = total_horizon_time - current_sim_time 

330 self.flex_model.do_step(t_start=0, t_sample=t_sample) 

331 

332 # save output 

333 for output in self.var_ref.outputs: 

334 self.flex_results.loc[ 

335 (self.env.now, current_sim_time + t_sample), output 

336 ] = self.flex_model.get_output(output).value 

337 

338 

339class FlexibilityShadowMINLPMPCConfig(minlp_mpc.MINLPMPCConfig): 

340 

341 baseline_input_names: list[str] = Field(default=[]) 

342 custom_input_names: list[Dict] = Field(default=[]) 

343 baseline_state_names: list[str] = Field(default=[]) 

344 full_control_names: list[str] = Field(default=[]) 

345 

346 baseline_agent_id: str = "" 

347 

348 casadi_sim_time_step: int = Field( 

349 default=0, 

350 description="Time step for simulation with Casadi simulator. " 

351 "Value is read from FlexQuantConfig", 

352 ) 

353 power_variable_name: str = Field( 

354 default=None, description="Name of the power variable in the " 

355 "shadow mpc model." 

356 ) 

357 storage_variable_name: Optional[str] = Field( 

358 default=None, description="Name of the storage variable in the " 

359 "shadow mpc model." 

360 ) 

361 

362 

363class FlexibilityShadowMINLPMPC(minlp_mpc.MINLPMPC): 

364 """Shadow MINLP-MPC for calculating positive/negatives flexibility offers. 

365 

366 """ 

367 

368 config: FlexibilityShadowMINLPMPCConfig 

369 

370 def __init__(self, *args, **kwargs): 

371 # initialize flex_results with None 

372 self.flex_results = None 

373 

374 super().__init__(*args, **kwargs) 

375 

376 # setup look up dict to track incoming inputs and states 

377 # (maps name as str to actual AgentVariable) 

378 input_names_list = [var["name"] for var in self.config.custom_input_names] 

379 self._track_base_comm_vars_dict: Dict[str, Union[AgentVariable, None]] = {} 

380 for comm_var in self.config.inputs + self.config.states: 

381 if (comm_var.name in self.config.full_control_names or 

382 comm_var.name + base_vars_to_communicate_suffix in 

383 self.config.baseline_input_names or 

384 comm_var.name + base_vars_to_communicate_suffix in 

385 self.config.baseline_state_names or 

386 comm_var.name in input_names_list): 

387 comm_var.value = None 

388 self._track_base_comm_vars_dict[comm_var.name] = comm_var.copy(deep=True) 

389 # set up necessary components if simulation is enabled 

390 if self.config.casadi_sim_time_step > 0: 

391 # generate a separate simulation model for integration to ensure 

392 # the model used in MPC optimization remains unaffected 

393 self.flex_model = type(self.model)(dt=self.config.casadi_sim_time_step) 

394 # generate the filename for the simulation results 

395 self.res_file_flex = self.config.optimization_backend["results_file"].replace( 

396 "_flex", "_sim_flex" 

397 ) 

398 # clear the casadi simulator result at the first time step if already exists 

399 try: 

400 os.remove(self.res_file_flex) 

401 except FileNotFoundError: 

402 pass 

403 

404 def register_callbacks(self): 

405 for control_var in self.config.controls + self.config.binary_controls: 

406 self.agent.data_broker.register_callback( 

407 name=control_var.name + full_trajectory_suffix, 

408 alias=control_var.name + full_trajectory_suffix, 

409 callback=self.calc_flex_callback, 

410 source=Source(agent_id=self.config.baseline_agent_id, module_id=None) 

411 ) 

412 for base_inputs in self.config.baseline_input_names: 

413 self.agent.data_broker.register_callback( 

414 name=base_inputs.removesuffix(base_vars_to_communicate_suffix), # update MPC variable 

415 alias=base_inputs, 

416 callback=self.calc_flex_callback, 

417 source=Source(agent_id=self.config.baseline_agent_id, module_id=None) 

418 ) 

419 for custom_inputs in self.config.custom_input_names: 

420 self.agent.data_broker.register_callback( 

421 name=custom_inputs["name"], 

422 alias=custom_inputs["alias"], 

423 callback=self.calc_flex_callback 

424 ) 

425 for base_states in self.config.baseline_state_names: 

426 self.agent.data_broker.register_callback( 

427 name=base_states.removesuffix(base_vars_to_communicate_suffix), # update MPC variable 

428 alias=base_states, 

429 callback=self.calc_flex_callback, 

430 source=Source(agent_id=self.config.baseline_agent_id, module_id=None) 

431 ) 

432 

433 super().register_callbacks() 

434 

435 def calc_flex_callback(self, inp: AgentVariable, name: str): 

436 """Ensure that all control trajectories and Baseline inputs/states 

437 have been set before starting the calculation. 

438 

439 """ 

440 # during provision do not calculate flex 

441 if self.get(glbs.PROVISION_VAR_NAME).value: 

442 return 

443 

444 # do not trigger callback on self set variables 

445 if self.agent.config.id == inp.source.agent_id: 

446 return 

447 

448 # get the value of the input 

449 vals = inp.value 

450 

451 if inp.name in self.config.full_control_names: 

452 if vals.isna().any(): 

453 vals = fill_nans(series=vals, method=MEAN) 

454 # set full controls to custom cia backend to constrain during market time 

455 if isinstance(self.optimization_backend, ConstrainedCasADiCIABackend): 

456 self.optimization_backend.config.full_controls_dict[inp.name] = vals 

457 # add time shift env.now to the mpc prediction index if it starts at t=0 

458 if vals.index[0] == 0: 

459 vals.index += self.env.time 

460 

461 # update value in the tracking dictionary 

462 self._track_base_comm_vars_dict[name].value = vals 

463 # set value 

464 self.set(name, vals) 

465 

466 # make sure all necessary inputs are set 

467 if all(x.value is not None for x in self._track_base_comm_vars_dict.values()): 

468 self.do_step() 

469 for _, comm_var in self._track_base_comm_vars_dict.items(): 

470 comm_var.value = None 

471 

472 def process(self): 

473 # the shadow mpc should only be run after the results of the baseline are sent 

474 yield self.env.event() 

475 

476 def set_output(self, solution): 

477 """Takes the solution from optimization backend and sends 

478 it to AgentVariables. 

479 

480 """ 

481 # Output must be defined in the config as "type"="pd.Series" 

482 if not self.config.set_outputs: 

483 return 

484 self.logger.info("Sending optimal output values to data_broker.") 

485 

486 # simulate with the casadi simulator 

487 self.sim_flex_model(solution) 

488 

489 df = solution.df 

490 if self.flex_results is not None: 

491 for output in self.var_ref.outputs: 

492 if output not in [ 

493 self.config.power_variable_name, 

494 self.config.storage_variable_name, 

495 ]: 

496 series = df.variable[output] 

497 self.set(output, series) 

498 # send the power and storage variable value from simulation results 

499 upsampled_output_power = self.flex_results[self.config.power_variable_name] 

500 self.set(self.config.power_variable_name, upsampled_output_power) 

501 if self.config.storage_variable_name is not None: 

502 upsampled_output_storage = self.flex_results[self.config.storage_variable_name] 

503 self.set(self.config.storage_variable_name, upsampled_output_storage.dropna()) 

504 else: 

505 for output in self.var_ref.outputs: 

506 series = df.variable[output] 

507 self.set(output, series) 

508 

509 def sim_flex_model(self, solution): 

510 """simulate the flex model over the preditcion horizon and save results 

511 

512 """ 

513 

514 # return if sim_time_step is not a positive integer and system is in provision 

515 if not (self.config.casadi_sim_time_step > 0 and not self.get(glbs.PROVISION_VAR_NAME).value): 

516 return 

517 

518 # read the defined simulation time step 

519 sim_time_step = self.config.casadi_sim_time_step 

520 mpc_time_step = self.config.time_step 

521 

522 # set the horizon length and the number of simulation steps 

523 total_horizon_time = int(self.config.prediction_horizon * self.config.time_step) 

524 n_simulation_steps = math.ceil(total_horizon_time / sim_time_step) 

525 

526 # read the current optimization result 

527 result_df = solution.df 

528 

529 # initialize the flex sim results Dataframe 

530 self._initialize_flex_results( 

531 n_simulation_steps, total_horizon_time, sim_time_step, result_df 

532 ) 

533 

534 # Update model parameters and initial states 

535 self._update_model_parameters() 

536 self._update_initial_states(result_df) 

537 

538 # Run simulation 

539 self._run_simulation( 

540 n_simulation_steps, sim_time_step, mpc_time_step, result_df, total_horizon_time 

541 ) 

542 

543 # set index of flex results to the same as mpc result 

544 store_results_df = self.flex_results.copy(deep=True) 

545 store_results_df.index = self.flex_results.index.tolist() 

546 

547 # save results 

548 if not os.path.exists(self.res_file_flex): 

549 store_results_df.to_csv(self.res_file_flex) 

550 else: 

551 store_results_df.to_csv(self.res_file_flex, mode="a", header=False) 

552 

553 # set the flex results format same as mpc result while updating Agentvariable 

554 self.flex_results.index = self.flex_results.index.get_level_values(1) 

555 

556 def _initialize_flex_results( 

557 self, n_simulation_steps, horizon_length, sim_time_step, result_df 

558 ): 

559 """Initialize the flex results dataframe with the correct dimension 

560 and index and fill with existing results from optimization 

561 

562 """ 

563 

564 # create MultiIndex for collocation points 

565 index_coll = pd.MultiIndex.from_arrays( 

566 [[self.env.now] * len(result_df.index), result_df.index], 

567 names=["time_step", "time"] 

568 # Match the names with multi_index but note they're reversed 

569 ) 

570 # create Multiindex for full simulation sample times 

571 index_full_sample = pd.MultiIndex.from_tuples( 

572 zip( 

573 [self.env.now] * (n_simulation_steps + 1), 

574 range(0, horizon_length + sim_time_step, sim_time_step), 

575 ), 

576 names=["time_step", "time"], 

577 ) 

578 # merge indexes 

579 new_index = index_coll.union(index_full_sample).sort_values() 

580 # initialize the flex results with correct dimension 

581 self.flex_results = pd.DataFrame(np.nan, index=new_index, columns=self.var_ref.outputs) 

582 

583 # Get the optimization outputs and create a series for fixed optimization outputs with the 

584 # correct MultiIndex format 

585 opti_outputs = result_df.variable[self.config.power_variable_name] 

586 fixed_opti_output = pd.Series( 

587 opti_outputs.values, 

588 index=index_coll, 

589 ) 

590 # fill the output value at the time step where it already exists in optimization output 

591 for idx in fixed_opti_output.index: 

592 if idx in self.flex_results.index: 

593 self.flex_results.loc[idx, self.config.power_variable_name] = fixed_opti_output[idx] 

594 

595 def _update_model_parameters(self): 

596 """update the value of module parameters with value from config, 

597 since creating a model just reads the value in the model class but not the config 

598 """ 

599 

600 for par in self.config.parameters: 

601 self.flex_model.set(par.name, par.value) 

602 

603 def _update_initial_states(self, result_df): 

604 """set the initial value of states""" 

605 

606 # get state values from the mpc optimization result 

607 state_values = result_df.variable[self.var_ref.states] 

608 # update state values with last measurement 

609 for state, value in zip(self.var_ref.states, state_values.iloc[0]): 

610 self.flex_model.set(state, value) 

611 

612 def _run_simulation( 

613 self, n_simulation_steps, sim_time_step, mpc_time_step, result_df, total_horizon_time 

614 ): 

615 """simulate with flex model over the prediction horizon""" 

616 

617 # get control and input values from the mpc optimization result 

618 control_values = result_df.variable[ 

619 [*self.var_ref.controls, *self.var_ref.binary_controls] 

620 ].dropna() 

621 input_values = result_df.parameter[self.var_ref.inputs].dropna() 

622 

623 # Get the simulation time step index 

624 sim_time_index = np.arange(0, (n_simulation_steps + 1) * sim_time_step, sim_time_step) 

625 

626 # Reindex the controls and inputs to sim_time_index 

627 control_values_full = control_values.copy().reindex(sim_time_index, method="ffill") 

628 input_values_full = input_values.copy().reindex(sim_time_index, method="nearest") 

629 

630 for i in range(0, n_simulation_steps): 

631 current_sim_time = i * sim_time_step 

632 

633 # Apply control and input values from the appropriate MPC step 

634 for control, value in zip( 

635 self.var_ref.controls, 

636 control_values_full.loc[current_sim_time, self.var_ref.controls], 

637 ): 

638 self.flex_model.set(control, value) 

639 

640 for binary_control, value in zip( 

641 self.var_ref.binary_controls, 

642 control_values_full.loc[current_sim_time, self.var_ref.binary_controls], 

643 ): 

644 self.flex_model.set(binary_control, value) 

645 

646 for input_var, value in zip( 

647 self.var_ref.inputs, input_values_full.loc[current_sim_time] 

648 ): 

649 # change the type of iterable input, since casadi model can't deal with iterable 

650 if issubclass(eval(self.flex_model.get(input_var).type), Iterable): 

651 self.flex_model.get(input_var).type = type(value).__name__ 

652 self.flex_model.set(input_var, value) 

653 

654 # do integration 

655 # reduce the simulation time step so that the total horizon time will not be exceeded 

656 if current_sim_time + sim_time_step <= total_horizon_time: 

657 t_sample = sim_time_step 

658 else: 

659 t_sample = total_horizon_time - current_sim_time 

660 self.flex_model.do_step(t_start=0, t_sample=t_sample) 

661 

662 # save output 

663 for output in self.var_ref.outputs: 

664 self.flex_results.loc[ 

665 (self.env.now, current_sim_time + t_sample), output 

666 ] = self.flex_model.get_output(output).value