Dynamic Operability Analysis of the HYPER Cyber-Physical Fuel Cell and Gas Turbine Hybrid System#
Author: Victor Alves - Carnegie Mellon University
HYPER (HYbrid PERformance) is a cyber-physical fuel cell and gas turbine hybrid power system at the National Energy Technology Laboratory (NETL). It converts coal into a high-purity carbon monoxide stream that powers a solid oxide fuel cell (SOFC) and can reduce the carbon footprint of coal-based power generation by up to \(90\%\) [11]. [11] uses it to demonstrate dynamic operability analysis of a linear time-invariant (LTI) system, illustrating:
the fast, exact dynamic achievable output funnel for an LTI system by state-space projection;
the dynamic Operability Index (dOI) at each time step, shown on the funnel;
robustness to disturbances via the intersection of achievable sets across scenarios;
an online funnel update for a perturbed initial state, giving the operability recovery horizon.
The identified HYPER model#
Transfer functions identified from the physical gas turbine-fuel cell system at NETL [11] relate the controlled outputs (turbine speed \(y_1\), cathode airflow \(y_2\)) to the manipulated inputs (electricity load \(u_1\), cold-air bypass \(u_2\)) and the disturbances (fuel valve \(w_1\), hot-air bypass \(w_2\)):
Nominal values [11]:
Variable |
Description |
Nominal |
|---|---|---|
\(y_1\) |
Turbine speed [rpm] |
40,500 |
\(y_2\) |
Cathode airflow [kg/min] |
45 (0.75 kg/s) |
\(u_1\) |
Electricity load [kW] |
40 |
\(u_2\) |
Cold-air bypass [%] |
40 |
\(w_1\) |
Fuel valve [%] |
50 |
\(w_2\) |
Hot-air bypass [%] |
10 |
[11] discretizes these transfer functions to a 4-state state-space model with a zero-order hold at \(\Delta t = 1\) min. The identified model is in deviation variables; the nominal references u_ref and y_ref keep all sets and plots in absolute units. [11] does not report the state-space matrices, so they are reconstructed here from the transfer functions, neglecting the sub-sample dead times (\(0.1\) to \(0.7\) min). (Cathode airflow is shown in kg/min, with nominal \(0.75\) kg/s \(= 45\) kg/min.)
Let’s start by importing the packages:
import time
import numpy as np
import matplotlib.pyplot as plt
# Plotly renderer that stays interactive in the compiled
# documentation website (loads plotly.js via require.js).
import plotly.io as pio
pio.renderers.default = 'notebook_connected'
Building the discrete-time matrices from the transfer functions: each first-order entry \(K/(s+a)\) contributes one state, discretized by zero-order hold as \(A = e^{-a \Delta t}\) and \(B = (1-e^{-a \Delta t})/a\).
Ts = 1.0 # Sampling time [min]
# Plant transfer functions: gains K and poles a, entry (i, j) = y_i / u_j
K = np.array([[-0.25, -0.065],
[-0.22, -1.43]])
a = np.array([[0.225, 0.125],
[0.69, 1.43]])
# Disturbance transfer functions, entry (i, j) = y_i / w_j
Kw = np.array([[0.17, 0.03],
[0.06, -1.23]])
aw = np.array([[0.125, 0.46],
[0.16, 1.43]])
def discretize_first_order_tfs(gains, poles):
"""
Discretize a 2x2 matrix of first-order transfer functions into a single
discrete-time state-space model.
Each entry of the matrix is a first-order transfer function
``gains[i, j] / (s + poles[i, j])``. It is realized with one state and
discretized with a zero-order hold at the sampling time ``Ts``, so the
full 2x2 matrix yields a four-state model ``x(k+1) = A x(k) + B u(k)``,
``y(k) = C x(k)``.
Parameters
----------
gains : numpy.ndarray, shape (2, 2)
Transfer-function gains; entry ``(i, j)`` is the gain from input ``j``
to output ``i``.
poles : numpy.ndarray, shape (2, 2)
Transfer-function poles (the ``a`` in ``K / (s + a)``), same layout.
Returns
-------
A : numpy.ndarray, shape (4, 4)
Discrete-time state matrix (``exp(-a * Ts)`` on the diagonal).
B : numpy.ndarray, shape (4, 2)
Discrete-time input matrix (``(1 - exp(-a * Ts)) / a`` per state).
C : numpy.ndarray, shape (2, 4)
Output matrix carrying the transfer-function gains.
"""
A = np.diag([np.exp(-poles[i, j] * Ts) for i in range(2) for j in range(2)])
B = np.zeros((4, 2))
C = np.zeros((2, 4))
s = 0
for i in range(2):
for j in range(2):
B[s, j] = (1 - np.exp(-poles[i, j] * Ts)) / poles[i, j]
C[i, s] = gains[i, j]
s += 1
return A, B, C
A, B, C = discretize_first_order_tfs(K, a) # 4-state plant model
Aw, Gw, Cw = discretize_first_order_tfs(Kw, aw) # 4-state disturbance model
# Sanity check: the discrete DC gain must match the continuous gains K/a
print('Discrete DC gain:')
print(np.round(C @ np.linalg.inv(np.eye(4) - A) @ B, 4))
print('Continuous DC gain (K/a):')
print(np.round(K / a, 4))
Discrete DC gain:
[[-1.1111 -0.52 ]
[-0.3188 -1. ]]
Continuous DC gain (K/a):
[[-1.1111 -0.52 ]
[-0.3188 -1. ]]
Nominal operating point and operability sets, in absolute units [11]: the AIS is electricity load \(30\) to \(50\) kW and cold-air bypass \(40\) to \(80\%\); the DOS is turbine speed \(40{,}490\) to \(40{,}510\) rpm and cathode airflow \(10\) to \(30\) kg/min.
# Nominal operating point (from the reference)
u_ref = np.array([40.0, 40.0]) # Electricity load [kW], cold-air bypass [%]
y_ref = np.array([40500.0, 45.0]) # Turbine speed [rpm], cathode airflow [kg/min]
w_ref = np.array([50.0, 10.0]) # Fuel valve [%], hot-air bypass [%]
# Available Input Set (AIS), absolute units
AIS_bounds = np.array([[30.0, 50.0], # Electricity load [kW]
[40.0, 80.0]]) # Cold-air bypass [%]
# Desired Output Set (DOS), absolute units
DOS_bounds = np.array([[40490.0, 40510.0], # Turbine speed [rpm]
[10.0, 30.0]]) # Cathode airflow [kg/min]
output_labels = ['Turbine speed [rpm]',
'Cathode airflow [kg/min]']
Mapping the nominal achievable output funnel with the fast linear method#
For an LTI system opyrability builds the funnel by state-space projection, propagating the reachable set forward exactly with no enumeration of input sequences [11]. Brute-force enumeration over 20 steps would need \(\sim 10^{12}\) simulations; the projection is exact and takes about a second [11].
Importing opyrability’s dynamic operability module:
from opyrability import dynamic_operability
Passing the model as an \(\{A, B, C\}\) dictionary triggers the fast linear method; u_ref and y_ref convert between the deviation-variable model and the absolute-unit sets. The dOI is evaluated against the DOS at each step and shown as the funnel coloring, and random (Monte Carlo) input sequences are overlaid as a check (each should fall inside the funnel). The figure is interactive: drag to rotate, scroll to zoom, right-drag to pan.
model = {'A': A, 'B': B, 'C': C, 'u_ref': u_ref, 'y_ref': y_ref}
start = time.time()
results = dynamic_operability(model,
np.zeros(4),
AIS_bounds,
DOS=DOS_bounds,
k_max=20,
monte_carlo=100,
seed=42,
engine='plotly',
labels=output_labels,
orientation='vertical')
print(f'Method selected automatically: {results["method"]}')
print(f'Elapsed time: {time.time() - start:.2f} s for 20 time steps')
print(f'dOI at each minute: {np.round(results["dOI"], 1)}')
Method selected automatically: state-space projection (linear)
Elapsed time: 1.12 s for 20 time steps
dOI at each minute: [ 29.4 63.3 77.3 83. 87.5 91.2 94.1 96.4 98.2 99.7 100. 100.
100. 100. 100. 100. 100. 100. 100. 100. ]
The funnel starts as the single nominal output point (dOI = 0) and grows as the inputs move the process: quickly along cathode airflow (fastest channel \(y_2/u_2\), \(1/1.43 \approx 0.7\) min) and slowly along turbine speed (slowest pole \(1/0.125 = 8\) min). Only \(29\%\) of the DOS is reachable at minute 1, and the full DOS at minute 11. Since the dOI is a property of the process, not of any controller, this limitation applies to every controller [11].
Assessing robustness to process disturbances#
The fuel valve and hot-air bypass disturbances are Gaussian about their nominal values with variances \(15\%\) of those values [11]. Here their \(95\%\) ranges (\(\pm 1.96\,\sigma\)) are gridded at the corners, each held as a sustained worst-case scenario; the robust outputs are the intersection of the scenario funnels. This corner treatment is a more conservative approximation of the ellipsoidal \(95\%\)-region analysis of [11]. The \(5\%\)-variance measurement noises of [11] are not included here.
Importing opyrability’s scenario evaluation module:
from opyrability import dynamic_operability_scenarios
The disturbance enters through an augmented 8-state model (4 plant states and 4 disturbance states). Each scenario fixes its disturbance, so we supply one step model per scenario (absolute inputs in, absolute outputs out) and build the funnel by simulation:
# Augmented model: plant states stacked with disturbance states
A8 = np.block([[A, np.zeros((4, 4))], [np.zeros((4, 4)), Aw]])
B8 = np.vstack([B, np.zeros((4, 2))])
G8 = np.vstack([np.zeros((4, 2)), Gw])
C8 = np.hstack([C, Cw])
# 95% disturbance ranges: variance = 15% of the nominal values (50 %, 10 %)
sigma_w = np.sqrt(0.15 * w_ref)
w_95 = 1.96 * sigma_w
print(f'95% disturbance ranges: fuel valve {w_ref[0]} +/- {w_95[0]:.2f} %, '
f'hot-air bypass {w_ref[1]} +/- {w_95[1]:.2f} %')
def make_step(w):
"""
Build the step model for a fixed disturbance ``w``.
Returns a function ``step(x, u) -> (x_next, y)`` that advances the
augmented (plant plus disturbance) state by one minute with the
disturbance held constant at ``w``. Inputs ``u`` and outputs ``y`` are in
absolute units; the underlying model is in deviation variables, so the
nominal references are subtracted from the inputs and added back to the
outputs.
Parameters
----------
w : array-like, shape (2,)
Sustained disturbance ``[fuel valve %, hot-air bypass %]``.
Returns
-------
callable
The step model ``step(x, u)``.
"""
w_dev = np.asarray(w, dtype=float) - w_ref
def step(x, u):
u_dev = np.asarray(u, dtype=float) - u_ref
x_next = A8 @ np.asarray(x, float) + B8 @ u_dev + G8 @ w_dev
return x_next, C8 @ x_next + y_ref
return step
def make_initial_state(w):
"""
Return the starting augmented state for disturbance ``w``.
The reactor starts at the nominal steady state, which is the zero vector
in deviation variables, independent of ``w``.
Parameters
----------
w : array-like, shape (2,)
Sustained disturbance (unused; present for a uniform interface).
Returns
-------
numpy.ndarray, shape (8,)
The initial augmented state.
"""
return np.zeros(8)
corner_scenarios = {
f'w = ({w_ref[0] + i * w_95[0]:.1f}, {w_ref[1] + j * w_95[1]:.1f}) %':
np.array([w_ref[0] + i * w_95[0], w_ref[1] + j * w_95[1]])
for i in (1, -1) for j in (1, -1)}
scenario_results = dynamic_operability_scenarios(make_step,
make_initial_state,
AIS_bounds,
scenarios=corner_scenarios,
DOS=DOS_bounds,
k_max=20,
method='nstep',
engine='plotly',
labels=output_labels,
orientation='vertical')
print('Robust dOI at each minute:',
np.round(scenario_results['intersection']['dOI'], 1))
95% disturbance ranges: fuel valve 50.0 +/- 5.37 %, hot-air bypass 10.0 +/- 2.40 %
Robust dOI at each minute: [ 7.3 26. 41.8 53.6 56. 57.4 58.4 59. 59.3 59.5 59.5 59.4 59.2 59. ]
Each sustained disturbance corner shifts the funnel (colored outlines); the blue region is their intersection, the outputs guaranteed for any disturbance realization. The guaranteed dOI plateaus at about \(60\%\) of the DOS. With disturbances varying in time inside the \(95\%\) ellipsoid, [11] recovers the entire DOS after \(8\) minutes; the corner treatment here is the conservative bound on that result.
Updating the funnel online for a perturbed initial state#
In operation a disturbance can push the process off the nominal state, so the funnel is rebuilt from the new state. Following [11], the initial state is moved to the steady state with outputs \(40{,}588\) rpm and \(133\) kg/min. The output funnel is a linear image of the state funnel, so the update is shown first in state space. The state variables have no physical meaning (they are realization-dependent), so the shapes differ from [11] while the behavior matches: HYPER is stable, so the updated funnel converges back to the nominal one, which sets the longest recovery horizon. The update needs no re-simulation; the stored funnel inequalities only shift their right-hand sides [11].
Importing opyrability’s state funnel visualization and online update modules:
from opyrability import plot_state_funnel, update_dynamic_funnel
Computing the perturbed initial state and comparing the nominal and perturbed state funnels (projected onto the first two state variables):
# Steady state consistent with outputs (40588 rpm, 133 kg/min)
y_perturbed = np.array([40588.0, 133.0])
dc_gain = C @ np.linalg.inv(np.eye(4) - A) @ B
u_ss = np.linalg.solve(dc_gain, y_perturbed - y_ref)
x0_perturbed = np.linalg.solve(np.eye(4) - A, B @ u_ss)
print(f'Perturbed initial outputs: {np.round(C @ x0_perturbed + y_ref, 1)}')
# State funnels from the nominal and the perturbed initial states
res_nominal = dynamic_operability(model,
np.zeros(4),
AIS_bounds,
k_max=20,
plot=False)
start = time.time()
res_perturbed = update_dynamic_funnel(res_nominal,
x0_perturbed)
t_update = time.time() - start
start = time.time()
res_recomputed = dynamic_operability(model,
x0_perturbed,
AIS_bounds,
k_max=20,
plot=False)
t_offline = time.time() - start
print(f'Online update (hyperplane shift): {1000 * t_update:.0f} ms, '
f'offline reconstruction: {t_offline:.2f} s')
funnels = {'Nominal initial state': res_nominal,
'Perturbed initial state': res_perturbed}
fig_states, _ = plot_state_funnel(funnels,
dims=(0, 1),
engine='plotly')
Perturbed initial outputs: [40588. 133.]
Online update (hyperplane shift): 2 ms, offline reconstruction: 1.00 s
The two state funnels start apart and merge over time. In the output space, the updated funnel starts at the perturbed point, away from the DOS, and travels toward the DOS box until it covers it:
results_updated = dynamic_operability(model,
x0_perturbed,
AIS_bounds,
DOS=DOS_bounds,
k_max=20,
monte_carlo=100,
seed=42,
engine='plotly',
labels=output_labels,
orientation='vertical')
print(f'dOI at each minute: {np.round(results_updated["dOI"], 1)}')
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
dOI at each minute: [ 0. 0. 0. 0. 0. 0. 23.4 50.5 73.4 92.8 100. 100.
100. 100. 100. 100. 100. 100. 100. 100. ]
The dOI is zero for the first \(6\) minutes: from the perturbed state, no desired output is reachable yet. The funnel reaches the DOS at minute \(7\) and covers it fully at minute \(11\). A complementary view intersects the updated funnel with the nominal one, showing when the nominally achievable outputs become reachable again:
def make_step_no_disturbance(scenario):
"""
Build the disturbance-free step model used for the initial-state
comparison.
The dynamics are identical for both scenarios (they differ only in their
starting state), so ``scenario`` does not change the returned model.
Parameters
----------
scenario : str
Scenario tag (unused; present for a uniform interface).
Returns
-------
callable
The step model ``step(x, u) -> (x_next, y)`` in absolute units.
"""
def step(x, u):
u_dev = np.asarray(u, dtype=float) - u_ref
x_next = A @ np.asarray(x, float) + B @ u_dev
return x_next, C @ x_next + y_ref
return step
def initial_state_for(scenario):
"""
Return the starting state for a scenario: the perturbed state for
``'perturbed'`` and the nominal steady state otherwise.
Parameters
----------
scenario : str
Either ``'perturbed'`` or ``'nominal'``.
Returns
-------
numpy.ndarray, shape (4,)
The initial state.
"""
return x0_perturbed.copy() if scenario == 'perturbed' else np.zeros(4)
init_scenarios = {'Nominal initial state': 'nominal',
'Perturbed initial state': 'perturbed'}
update_results = dynamic_operability_scenarios(make_step_no_disturbance,
initial_state_for,
AIS_bounds,
scenarios=init_scenarios,
DOS=DOS_bounds,
k_max=20,
method='nstep',
engine='plotly',
labels=output_labels)
overlap_update = update_results['intersection']['volumes']
first_k = int(np.argmax(overlap_update > 0)) + 1
print(f'First overlap between nominal and perturbed funnels: k = {first_k} min')
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
First overlap between nominal and perturbed funnels: k = 8 min
The intersection is empty for the first \(7\) minutes; the first overlap appears at minute \(8\), the recovery horizon of [11]. [11] notes this horizon can serve as the prediction horizon for model predictive control.
The with-disturbance funnel: Gaussian uncertainty via hyperplane shrinkage#
The corner treatment holds disturbances at worst-case values, which is conservative. [11] instead treats them as Gaussian sequences and shrinks each funnel face inward so the remaining outputs stay achievable for every realization in the \(95\%\) region. opyrability applies the single-stage output-space form: it propagates the disturbance covariance to an output covariance at each step and shrinks the updated funnel accordingly. The measurement noises of [11] are not included, as their reported variances cannot be reconciled with the rpm output scale.
Importing opyrability’s Gaussian robustness modules:
from opyrability import (propagate_output_covariance,
gaussian_robust_funnel, plot_dynamic_funnel)
Propagating the disturbance covariance \(\Sigma_w\) (\(15\%\)-variance) through the disturbance model and shrinking the updated funnel:
Sigma_w = np.diag(0.15 * w_ref)
Sigma_y_k = propagate_output_covariance(Aw,
Gw,
Cw,
Sigma_w,
k_max=20)
robust_updated = gaussian_robust_funnel(results_updated,
Sigma_y_k,
confidence=0.95,
DOS=DOS_bounds)
fig_rob, _ = plot_dynamic_funnel(robust_updated,
DOS=DOS_bounds,
dOI=robust_updated['dOI'],
labels=output_labels,
engine='plotly')
print('Robust dOI (Gaussian, updated funnel): ',
np.round(robust_updated['dOI'], 1))
print('Robust dOI (corner scenarios): ',
np.round(scenario_results['intersection']['dOI'], 1))
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
No intersection found between AS and DS. OI = 0.
Robust dOI (Gaussian, updated funnel): [ 0. 0. 0. 0. 0. 0. 13. 39.9 62.7 81.9 98.2 100.
100. 100. 100. 100. 100. 100. 100. 100. ]
Robust dOI (corner scenarios): [ 7.3 26. 41.8 53.6 56. 57.4 58.4 59. 59.3 59.5 59.5 59.4 59.2 59. ]
The Gaussian disturbances delay full operability by one minute: the entire DOS is guaranteed at minute \(12\), versus minute \(11\) without disturbances. The corner intersection never exceeds about \(60\%\): holding disturbances at worst-case values forever is far stronger than letting them fluctuate inside the \(95\%\) region.
Conclusions#
Key results:
The 20-minute HYPER funnel was computed in about a second by exact state-space projection, which
dynamic_operabilityselects automatically for an \(\{A, B, C\}\) model, versus \(\sim 10^{12}\) simulations for brute-force enumeration [11].The deviation-variable identified model was analyzed directly in absolute units via
u_refandy_ref.The dOI shows HYPER’s transient limit: the full DOS is reachable only at minute \(11\), for any controller [11].
Sustained worst-case disturbances keep about \(60\%\) of the DOS guaranteed; under Gaussian disturbances the full DOS returns at minute \(12\), one minute after the disturbance-free case.
The funnel is updated for a new initial state by an exact hyperplane shift [11] in milliseconds (versus about a second to recompute); the perturbed funnel reaches the DOS at minute \(7\), covers it at minute \(11\), and overlaps the nominal funnel at the recovery horizon of minute \(8\) [11].
Workflow:
Reconstruct the discrete-time LTI model from the transfer functions.
Set the AIS and DOS in absolute units with
u_refandy_ref.dynamic_operabilitymaps the funnel and dOI by the fast linear method.dynamic_operability_scenariosintersects funnels over disturbance and initial-state scenarios.plot_state_funnelcompares the nominal and updated state funnels.update_dynamic_funnelupdates online, andgaussian_robust_funnelshrinks the funnel for Gaussian uncertainty.