Concurrent simulation framework to enable feedback in system-level control

Concurrent simulation framework to enable feedback in system-level control

Change the H2Integrate framework so that the system-level model can be run in both sequential-type simulations and concurrent-type simulations. A sequential-type simulation is where the inputs and outputs are computed for the whole simulation period (8760 steps) on a subsystem-by subsystem basis:

A concurrent-type simulation is where the inputs and outputs are computed for all subsystems on a step-by-step basis:

Currently, only the sequential-type simulation is possible in H2Integrate, which makes it impractical to include feedback in system-level control decisions.

Proposed solution

This is a running code mock up of a simulation paradigm that could be used to run both sequential-type simulations and concurrent-type simulations in H2Integrate. This example relies on the definition of N_sim, the total number of steps in the simulation and N_step, the number of steps simulated in each computation step, both of which are made available throughout the model.

To run a sequential-type simulation set N_sim = 8760 and N_step = 8760 and call the openmdao prob.run_model(). Openmdao will run the model, computing the 8760-length inputs and outputs on a subsystem-by-subsystem basis as it traverses the system graph.

To run a concurrent-type simulation set N_sim = 8760 and N_step = 1 and call prob.run_model() in a loop while updating the time step index (see example). Openmdao will run the model, computing the 1-length inputs and outputs for each subsystem. When called in a loop, all subsystem models inputs and outputs will be computed concurrently for each step of the simulation.

import numpy as np
import matplotlib.pyplot as plt

import openmdao.api as om


class ModelBaseClass(om.ExplicitComponent):

    def initialize(self):

        # Number of steps in overall simulation, often 8760.
        self.options.declare(
            "N_sim", default=1, types=int, desc="number of steps in simulation"
        )

        # Number of steps to simulate for each compute call. This will be 8760 if running an annual
        # simulation and 1 if running a feedback steppable simulation.
        self.options.declare(
            "N_step", default=1, types=int, desc="number of steps per compute call"
        )

    def setup(self):

        # Add time step index as an input to all models
        self.add_discrete_input("k", val=0, desc="Time step index")


class Resource(ModelBaseClass):
    """
    A mock up of a resource class.
    """

    def setup(self):

        # Add time step index as input
        super().setup()

        self.add_input("mean_resource", shape=1)
        self.add_output("resource_signal", shape=(self.options["N_sim"],))

        self.use_cached_resource = False

    def get_resource_signal(self, mean_resource):
        """
        Create a synthetic resource signal with the mean value provided using a
        sine wave.
        """

        if self.use_cached_resource:

            return self.resource_signal

        else:
            # Generate synthetic resource signal. This is where the API call to the resource
            # database would happen.
            t = np.arange(0, self.options["N_sim"], 1)
            resource = mean_resource * np.ones(self.options["N_sim"])
            resource += (
                0.3 * mean_resource * np.sin(2 * np.pi / self.options["N_sim"] * t)
            )

            # Flip use cached flag to true so the resource signal is not re-generated every step
            self.use_cached_resource = True
            self.resource_signal = resource

            return resource

    def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):

        # Extract the current time step k and time step range for which this
        # compute call is relevant.
        k = discrete_inputs["k"]
        t_index = np.arange(k, k + self.options["N_step"], 1)

        # Compute/retreive or access the cached resource signal
        mean_resource = inputs["mean_resource"]
        resource_signal = self.get_resource_signal(mean_resource)

        # Send the resource signal for the relevant time step range to downstream systems.
        outputs["resource_signal"][t_index] = resource_signal[t_index]


class Converter(ModelBaseClass):
    """
    A mock up of a converter class
    """

    def setup(self):

        # add time step index as input
        super().setup()

        # Inputs and outputs are instantiated/stored as N_sim-length vectors and are filled
        # depending on the number of steps in the simulation. If it is an annual-type simulation,
        # then the whole 8760-length vector will be filled in one compute call. If it is a
        # steppable-type simulation then inputs and outputs will be filled progressively with N_step
        # length chunks of data every time compute is called.
        self.add_input("commodity_in", shape=(self.options["N_sim"],))
        self.add_output("commodity_out", shape=(self.options["N_sim"],))

    def compute_conversion(self, commodity_in):

        # Stand-in for a more complicated conversion calculation or the subsystem model contained in
        # the conversion subsystem wrapper.
        commodity_out = 0.5 * commodity_in

        return commodity_out

    def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):

        # Extract the current time step k and time step range for which this
        # compute call is relevant.
        k = discrete_inputs["k"]
        t_index = np.arange(k, k + self.options["N_step"], 1)

        # Extract subsystem inputs for the relevant time period.
        commodity_in = inputs["commodity_in"][t_index]

        # Run the model contained in the openmdao wrapper.
        commodity_out = self.compute_conversion(commodity_in)

        # Fill the output vector for the relevant time range.
        outputs["commodity_out"][t_index] = commodity_out


class Demand(ModelBaseClass):
    """
    A mock up of a resource class.
    """

    def setup(self):

        # add time step index as input
        super().setup()

        self.add_input("commodity_in", shape=(self.options["N_sim"],))
        self.add_input("demand_signal", shape=(self.options["N_sim"],))
        self.add_output("residual", shape=(self.options["N_sim"],))

    def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):

        # Extract the current time step k and time step range for which this
        # compute call is relevant.
        k = discrete_inputs["k"]
        t_index = np.arange(k, k + self.options["N_step"], 1)

        outputs["residual"][t_index] = (
            inputs["commodity_in"][t_index] - inputs["demand_signal"][t_index]
        )


class SystemAnnual(om.Group):
    """
    Group organization/construction for annual-type simulation run
    """

    def setup(self):

        # Simulation runs for 10 steps.
        N_sim = 10

        # All 10 steps are computed at once.
        N_step = 10

        self.add_subsystem(
            "resource", Resource(N_sim=N_sim, N_step=N_step), promotes_inputs=["k"]
        )
        self.add_subsystem(
            "converter", Converter(N_sim=N_sim, N_step=N_step), promotes_inputs=["k"]
        )
        self.add_subsystem(
            "demand", Demand(N_sim=N_sim, N_step=N_step), promotes_inputs=["k"]
        )

        self.connect("resource.resource_signal", "converter.commodity_in")
        self.connect("converter.commodity_out", "demand.commodity_in")


class SystemSteppable(om.Group):
    """
    Group organization/construction for annual-type simulation run
    """

    def setup(self):

        # Simulation runs for 10 steps.
        N_sim = 10

        # Only 1 step is computed at a time.
        N_step = 1

        self.add_subsystem(
            "resource", Resource(N_sim=N_sim, N_step=N_step), promotes_inputs=["k"]
        )
        self.add_subsystem(
            "converter", Converter(N_sim=N_sim, N_step=N_step), promotes_inputs=["k"]
        )
        self.add_subsystem(
            "demand", Demand(N_sim=N_sim, N_step=N_step), promotes_inputs=["k"]
        )

        self.connect("resource.resource_signal", "converter.commodity_in")
        self.connect("converter.commodity_out", "demand.commodity_in")


# Skip one or the other for troubleshooting
run_flag = dict(annual=True, steppable=True)


if run_flag["annual"]:

    # build the annual model and problem
    prob_annual = om.Problem()
    sys_annual = SystemAnnual()

    prob_annual.model.add_subsystem("sys", sys_annual)

    prob_annual.setup()

    prob_annual.set_val("sys.resource.mean_resource", 20)
    prob_annual.set_val("sys.demand.demand_signal", 10 * np.ones(10))

    # Run the model in one shot
    prob_annual.run_model()

    inputs_annual = prob_annual.model.list_inputs()
    outputs_annual = prob_annual.model.list_outputs()


if run_flag["steppable"]:

    # build the steppable model and problem
    prob_steppable = om.Problem()
    sys_steppable = SystemSteppable()

    prob_steppable.model.add_subsystem("sys", sys_steppable)

    prob_steppable.setup()

    prob_steppable.set_val("sys.resource.mean_resource", 20)
    prob_steppable.set_val("sys.demand.demand_signal", 10 * np.ones(10))

    # Run the model in a loop
    for k in range(10):

        # Update time index
        prob_steppable.set_val("sys.k", k)

        # Run the next step
        prob_steppable.run_model()

    """ Note
    I'm sure there is a better location for this loop than in the user-level run script. Maybe in 
    `H2IntegrateModel.run()`?
    """

    inputs_steppable = prob_steppable.model.list_inputs()
    outputs_steppable = prob_steppable.model.list_outputs()


n_puts = max(len(inputs_annual), len(outputs_annual))
# Plot inputs and outputs from simulation for comparison and sanity check
fig, ax = plt.subplots(n_puts, 2, sharex="all", layout="constrained", figsize=(10, 8))


def plot_input_output(ioputs, colnum, order):

    for i in range(len(order)):
        for name, val in ioputs:
            if name == order[i]:
                break

        iop = val["val"]

        ax[i, colnum].plot(iop)
        ax[i, colnum].set_title(name)


inp_order = [io[0] for io in inputs_annual]
outp_order = [io[0] for io in outputs_annual]

plot_input_output(inputs_annual, 0, order=inp_order)
plot_input_output(outputs_annual, 1, order=outp_order)

plot_input_output(inputs_steppable, 0, order=inp_order)
plot_input_output(outputs_steppable, 1, order=outp_order)

Alternatives considered

Additional context

Related to #204

Start by reviewing the current simulation execution path and the possible H2IntegrateModel.run() entry point mentioned in the issue. Compare it with the provided annual and steppable OpenMDAO mock-up, then identify how both modes and feedback would be represented. Done means H2Integrate supports sequential and concurrent simulations with matching outputs and a usable feedback path.