NMPC for rSOC example

README-developer.md

@@ -23,7 +23,7 @@ Clone the repository from GitHub, set up your Python environment as you usually
 pip install -r requirements-dev.txt
 ```
 
-Note: if you have IDAES installed in your current environment, it will uninstall it and install the latest version from the main branch on Github. You can run `pip uninstall idaes` and reinstall it from your local repository if you need to test examples against a local branch of IDAES.
+Note: if you have IDAES installed in your current environment, it will uninstall it and install the latest version from the main branch on Github. You can run `pip uninstall idaes-pse` and reinstall it from your local repository if you need to test examples against a local branch of IDAES.
 
 The configuration of the installation is stored in `pyproject.toml`.
 

idaes_examples/mod/power_gen/soc_nmpc/nmpc_tracking_LMP_settings.py

@@ -0,0 +1,121 @@
+import numpy as np
+import pandas as pd
+
+import pyomo.environ as pyo
+
+import idaes
+from idaes.core.solvers import get_solver
+from idaes.core.solvers import use_idaes_solver_configuration_defaults
+
+# Set up time discretization.
+pred_step_num = 5
+t_step = 75.0  # s
+pred_horizon = pred_step_num * t_step
+
+t_start = 1 * 30 * 60
+t_end = 2 * 60 * 60
+
+df_full = pd.read_csv("operations_new.csv", index_col=0)
+
+df = df_full
+df.reset_index(inplace=True)
+df.loc[:, 'hour'] = df.index
+df.set_index('hour', inplace=True)
+
+n_hour_point, n_vars = df.shape
+t_ramp = 10 * 60
+t_settle = 50 * 60
+
+dt_set = [t_start]
+scenario_names = ['t_start']
+for k in range(n_hour_point - 1):
+    dt_set.append(t_ramp)
+    dt_set.append(t_settle)
+    scenario_names.append('t_ramp')
+    scenario_names.append('t_settle')
+
+time_mark_list = [sum(dt_set[:j]) for j in range(len(dt_set) + 1)]
+sim_horizon = time_mark_list[-1] #- (t_settle - t_end)
+
+sim_nfe = int((sim_horizon) / t_step)
+sim_time_set = np.linspace(0, sim_nfe*t_step, sim_nfe+1)
+traj_time_set = np.linspace(0, sim_nfe*t_step+pred_horizon, sim_nfe+pred_step_num+1)
+# traj_time_set = np.linspace(0, sim_horizon+pred_horizon, int((sim_horizon+pred_horizon)/50)+1)
+
+# %% Set up ipopt.
+use_idaes_solver_configuration_defaults()
+idaes.cfg.ipopt.options.nlp_scaling_method = "gradient-based"
+idaes.cfg.ipopt.options.halt_on_ampl_error = "no"
+#linear_scaling_on_demand = no
+# idaes.cfg.ipopt_l1.options.OF_restoration_method = "l1"
+# idaes.cfg.ipopt_l1.options.OF_l1_epsilon = 0.5
+idaes.cfg.ipopt.options.max_iter = 1000
+idaes.cfg.ipopt.options.linear_solver = "ma57"
+idaes.cfg.ipopt.options.OF_ma57_pivtol = 1e-06  # default = 1e-08
+idaes.cfg.ipopt.options.OF_ma57_automatic_scaling = "yes"
+idaes.cfg.ipopt.options.OF_print_info_string = "yes"
+idaes.cfg.ipopt.options.OF_warm_start_init_point = "yes"
+idaes.cfg.ipopt.options.OF_warm_start_mult_bound_push = 1e-06
+idaes.cfg.ipopt.options.OF_warm_start_bound_push = 1e-06
+idaes.cfg.ipopt.options.tol = 1e-04  # default = 1e-08
+idaes.cfg.ipopt.options.mu_init = 1e-05  # default = 1e-01
+idaes.cfg.ipopt.options.bound_push = 1e-06  # default = 1e-02
+idaes.cfg.ipopt.options.bound_relax_factor = 1e-06  # default = 1e-08
+
+# idaes.cfg.ipopt.options.constr_viol_tol = 1e-3 # default = 1e-4
+solver_v1 = pyo.SolverFactory('ipopt')
+
+solver_options = {
+    "nlp_scaling_method": "gradient-based",
+    "halt_on_ampl_error": "no",
+    "max_iter": 1000,
+    "linear_solver": "ma57",
+    "ma57_pivtol": 1e-6,
+    "ma57_automatic_scaling": "yes",
+    "print_info_string": "no",
+    "warm_start_init_point": "yes",
+    "warm_start_mult_bound_push": 1e-6,
+    "warm_start_bound_push": 1e-6,
+    "tol": 1e-4,
+    "mu_init": 1e-5,
+    "bound_push": 1e-6,
+    "bound_relax_factor": 1e-6
+}
+solver_v2 = get_solver("ipopt_v2", options=solver_options, writer_config={"linear_presolve": False, "scale_model": True})
+
+
+def pass_all_vars_plant(target_model, source_model):
+    def get_all_vars(m):
+        all_vars = []
+        for var in m.component_objects(pyo.Var):
+            if isinstance(var, pyo.Var):
+                if var.name not in ('fs.p1', 'fs.n1',
+                                    'fs.p2', 'fs.n2',
+                                    'fs.p3', 'fs.n3',
+                                    'fs.soc_module.solid_oxide_cell.fuel_electrode.p',
+                                    'fs.soc_module.solid_oxide_cell.fuel_electrode.n',
+                                    'fs.p', 'fs.n',
+                                    'fs.avg_curr_density_p', 'fs.avg_curr_density_n',
+                                    'fs.pointwise_curr_density_p', 'fs.pointwise_curr_density_n',
+                                    # var in the controller but not in the plant
+                                    ):
+                    all_vars.append(var)
+        all_vars_temp = [v for v in all_vars if m.fs.time in v.index_set().subsets()]
+        return all_vars_temp
+
+    def set_var_ics(var_target, var_source):
+        for idxs, v in var_target.items():
+            if v.fixed:
+                continue
+            if type(idxs) == float:
+                v.set_value(var_source[idxs].value)
+            elif type(idxs) == tuple:
+                v.set_value(var_source[tuple([idxs[0], *idxs[1:]])].value)
+        return None
+
+    target_model_vars = get_all_vars(target_model)
+    source_model_vars = get_all_vars(source_model)
+    for var_target, var_source in zip(target_model_vars, source_model_vars):
+        set_var_ics(var_target, var_source)
+
+

idaes_examples/mod/power_gen/soc_nmpc/soec_attributes.py

@@ -0,0 +1,166 @@
+import numpy as np
+
+import pyomo.environ as pyo
+from pyomo.core.base.componentuid import ComponentUID
+
+from idaes_examples.mod.power_gen.soc_nmpc.nmpc_tracking_LMP_settings import time_mark_list, scenario_names, df
+
+class SocAttr:
+    def __init__(self):
+        self.alias_dict = {
+            "fs.lmp": "lmp",
+            "fs.total_electric_power": "total_electric_power",
+            "fs.h2_mass_production": "h2_production_rate",
+            "fs.h2_mass_consumption": "h2_consumption_rate",
+            "fs.soc_module.potential_cell": "potential",
+            "fs.soc_module.fuel_outlet_mole_frac_comp_H2": "soc_fuel_outlet_mole_frac_comp_H2",
+            "fs.makeup_mix._flow_mol_makeup_ref": "makeup_feed_rate",
+            "fs.sweep_blower._flow_mol_inlet_ref": "sweep_feed_rate",
+            "fs.feed_heater.electric_heat_duty": "feed_heater_duty",
+            "fs.feed_heater._temperature_outlet_ref": "feed_heater_outlet_temperature",
+            "fs.soc_module._temperature_fuel_outlet_ref": "fuel_outlet_temperature",
+            "fs.sweep_heater.electric_heat_duty": "sweep_heater_duty",
+            "fs.sweep_heater._temperature_outlet_ref": "sweep_heater_outlet_temperature",
+            "fs.soc_module._temperature_oxygen_outlet_ref": "sweep_outlet_temperature",
+            "fs.stack_core_temperature": "stack_core_temperature",
+            "fs.feed_recycle_split.recycle_ratio": "fuel_recycle_ratio",
+            "fs.sweep_recycle_split.recycle_ratio": "sweep_recycle_ratio",
+            "fs.sweep_recycle_split.mixed_state.mole_frac_comp[O2]": "oxygen_out",
+            "fs.feed_recycle_mix.mixed_state.mole_frac_comp[H2]": "hydrogen_in",
+            "fs.condenser_split.recycle_ratio": "vgr_recycle_ratio",
+            "fs.makeup_mix.makeup_mole_frac_comp_H2": "makeup_mole_frac_comp_H2",
+        }
+
+        self.setpoint_params = {
+            "fs.lmp":"lmp",
+            "fs.total_electric_power_sp": "total_electric_power",
+            "fs.h2_target_sp": "h2_production_rate",
+            "fs.h2_consumption_target_sp": "h2_consumption_rate",
+            "fs.potential_sp": "potential",
+            "fs.soc_fuel_outlet_mole_frac_comp_H2_sp": "soc_fuel_outlet_mole_frac_comp_H2",
+            "fs.makeup_feed_rate_sp": "makeup_feed_rate",
+            "fs.sweep_feed_rate_sp": "sweep_feed_rate",
+            "fs.feed_heater_duty_sp": "feed_heater_duty",
+            "fs.feed_heater_outlet_temperature_sp": "feed_heater_outlet_temperature",
+            "fs.fuel_outlet_temperature_sp": "fuel_outlet_temperature",
+            "fs.sweep_heater_duty_sp": "sweep_heater_duty",
+            "fs.sweep_heater_outlet_temperature_sp": "sweep_heater_outlet_temperature",
+            "fs.sweep_outlet_temperature_sp": "sweep_outlet_temperature",
+            "fs.stack_core_temperature_sp": "stack_core_temperature",
+            "fs.fuel_recycle_ratio_sp": "fuel_recycle_ratio",
+            "fs.sweep_recycle_ratio_sp": "sweep_recycle_ratio",
+            "fs.oxygen_out_sp": "oxygen_out",
+            "fs.hydrogen_in_sp": "hydrogen_in",
+            "fs.vgr_recycle_ratio_sp": "vgr_recycle_ratio",
+            "fs.makeup_mole_frac_comp_H2_sp": "makeup_mole_frac_comp_H2",
+        }
+        return None
+
+    def make_sp_trajs(self, tset):
+        """
+        Function to create setpoint trajectories for all variables
+        """
+        def make_one_var_traj(alias, tset):
+            """
+            Function to write a trajectory for a variable 
+            """
+            def make_traj_point(t, t1, t2, y1, y2):
+                """
+                Interpolation between two time points
+                """
+                a = (y2 - y1) / (t2 - t1)
+                return y1 + (t - t1) * a
+
+            var_target = {t: 0.0 for t in tset}
+            for t in tset:
+                for i in range(len(time_mark_list)):
+                    if t < time_mark_list[i]:
+                        lmp_sp_index = np.floor(i / 2) 
+                        scenario_index = i - 1
+                        current_scenario = scenario_names[scenario_index]
+
+                        if current_scenario == "t_ramp":
+                            var_target[t] = make_traj_point(t,
+                                                            time_mark_list[i - 1], time_mark_list[i],
+                                                            df[alias][lmp_sp_index - 1], df[alias][lmp_sp_index])
+                        elif current_scenario == "t_settle":
+                            var_target[t] = df[alias][lmp_sp_index]
+                        elif current_scenario == "t_start":
+                            var_target[t] = df[alias][0]
+                        break
+                    else:
+                        var_target[t] = df[alias].iloc[-1]
+
+            return var_target
+
+        var_trajs = {alias: {} for alias in self.alias_dict.values()}
+        for alias in var_trajs:
+            var_trajs[alias] = make_one_var_traj(alias, tset)
+
+        return var_trajs
+
+    def update_sp(self, m , t_base, var_targets):
+        for key, val in self.setpoint_params.items():
+            component = ComponentUID(key).find_component_on(m)
+            component_targets = var_targets[val]
+            for t in m.fs.time:
+                component[t].value = component_targets[t_base + t]
+
+    def update_coeff(self, m):
+        def adapt_sp_order(var):
+            if np.sqrt(pyo.value(var)) <= 1e-5:
+                return 1e+5
+            nearest_order = int(np.round(np.log10(np.sqrt(pyo.value(var)))))
+            return np.round(np.power(1e-01, nearest_order), decimals=nearest_order)
+
+        for t in m.fs.time:
+            m.fs.makeup_feed_rate_coeff[t] = adapt_sp_order(
+                (m.fs.makeup_mix.makeup.flow_mol[t] - m.fs.makeup_feed_rate_sp[t]) ** 2
+            )
+            m.fs.sweep_feed_rate_coeff[t] = adapt_sp_order(
+                (m.fs.sweep_blower.inlet.flow_mol[t] - m.fs.sweep_feed_rate_sp[t]) ** 2
+            )
+            m.fs.potential_coeff[t] = adapt_sp_order(
+                (m.fs.soc_module.potential_cell[t] - m.fs.potential_sp[t]) ** 2
+            )
+            m.fs.fuel_recycle_ratio_coeff[t] = adapt_sp_order(
+                (m.fs.feed_recycle_split.recycle_ratio[t] - m.fs.fuel_recycle_ratio_sp[t]) ** 2
+            )
+            m.fs.sweep_recycle_ratio_coeff[t] = adapt_sp_order(
+                (m.fs.sweep_recycle_split.recycle_ratio[t] - m.fs.sweep_recycle_ratio_sp[t]) ** 2
+            )
+            m.fs.vgr_recycle_ratio_coeff[t] = adapt_sp_order(
+                (m.fs.condenser_split.recycle_ratio[t] - m.fs.vgr_recycle_ratio_sp[t]) ** 2
+            )
+
+            m.fs.makeup_mole_frac_comp_H2_coeff[t] = adapt_sp_order(
+                (m.fs.makeup_mix.makeup_mole_frac_comp_H2[t] - m.fs.makeup_mole_frac_comp_H2_sp[t]) ** 2
+            )
+
+            m.fs.soc_fuel_outlet_mole_frac_comp_H2_coeff[t] = adapt_sp_order(
+                (m.fs.soc_module.fuel_outlet_mole_frac_comp_H2[t] - m.fs.soc_fuel_outlet_mole_frac_comp_H2_sp[t]) ** 2
+            )
+            m.fs.feed_heater_outlet_temperature_coeff[t] = adapt_sp_order(
+                (m.fs.feed_heater.outlet.temperature[t] - m.fs.feed_heater_outlet_temperature_sp[t]) ** 2
+            )
+            m.fs.sweep_heater_outlet_temperature_coeff[t] = adapt_sp_order(
+                (m.fs.sweep_heater.outlet.temperature[t] - m.fs.sweep_heater_outlet_temperature_sp[t]) ** 2
+            )
+            m.fs.fuel_outlet_temperature_coeff[t] = adapt_sp_order(
+                (m.fs.soc_module.fuel_outlet.temperature[t] - m.fs.fuel_outlet_temperature_sp[t]) ** 2
+            )
+            m.fs.sweep_outlet_temperature_coeff[t] = adapt_sp_order(
+                (m.fs.soc_module.oxygen_outlet.temperature[t] - m.fs.sweep_outlet_temperature_sp[t]) ** 2
+            )
+            m.fs.stack_core_temperature_coeff[t] = adapt_sp_order(
+                (m.fs.stack_core_temperature[t] - m.fs.stack_core_temperature_sp[t]) ** 2
+            )
+            m.fs.oxygen_out_coeff[t] = adapt_sp_order(
+                (m.fs.sweep_recycle_split.mixed_state[t].mole_frac_comp['O2'] - m.fs.oxygen_out_sp[t]) ** 2
+            )
+            m.fs.hydrogen_in_coeff[t] = adapt_sp_order(
+                (m.fs.feed_recycle_mix.mixed_state[t].mole_frac_comp['H2'] - m.fs.hydrogen_in_sp[t]) ** 2
+            )
+
+
+