Tutorial 6: State-space, Kalman and LPV export
Goal: turn the linearizations of a sweep into the models an estimator or a controller needs: the state-space model of one operating point, a reduced model of its structural modes, a Kalman filter that reconstructs a load from two accelerometers, a mode-tracked family scheduled on wind speed, and a family of named physical states.
Script: examples/06_state_space_export.py · Data: a synthetic sweep of ten
operating points, 4 to 15 m/s (see The synthetic turbine) ·
Run time: seconds
The code
from vane import ModalPipeline, discover_operating_points
from vane.export import (
build_physical_lpv_family,
load_lpv_family,
save_lpv_family,
write_rom,
)
from vane.sysid import (
CovarianceSet,
StateSpace,
build_kalman_model,
build_rom_family,
discretize_process_noise,
filter_consistency,
load_rom_family,
measurement_noise_from_specs,
modal_process_noise,
modal_reduction,
reconstruct_load,
run_filter,
save_rom_family,
schedule_model,
)
WIND_SPEEDS = (4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 11.0, 13.0, 15.0) # m/s
BAND_HZ = (0.1, 2.0)
SENSORS = ("YawBrTAxp", "YawBrTAyp")
SENSOR_STD = (0.02, 0.02) # m/s^2
DT = 0.05 # s
N_STEPS = 6000 # 300 s
def main(argv: list[str] | None = None) -> int:
"""Run the tutorial.
Parameters
----------
argv : list[str] or None, optional
Command-line arguments; ``None`` reads ``sys.argv``.
Returns
-------
int
The exit status.
"""
args = output_parser(__doc__).parse_args(argv)
with output_directories(args.output) as (data, figures):
# 1. The sweep and the modal pipeline on the wind-speed axis.
SyntheticTurbine().write_sweep(data / "sweep", WIND_SPEEDS, n_azimuths=12)
points = discover_operating_points(data / "sweep")
result = ModalPipeline().run(points, parameter_name="wind_speed")
winds = [mbc.wind_speed for mbc in result.mbc_results]
# 2. The full-order model at 8 m/s.
index = winds.index(8.0)
system = result.state_space(index)
print(
f"State space at {winds[index]} m/s: {system.n_states} states, "
f"{system.n_inputs} inputs, {system.n_outputs} outputs, "
f"stable: {system.is_stable()}"
)
print(f" inputs: {system.input_names}")
print(f" outputs: {system.output_names}")
# 3. Modal reduction to the structural modes that drive the tower-base load.
rom = modal_reduction(
system, n_modes=4, frequency_band_hz=BAND_HZ, rank_output="TwrBsMyt"
)
reduced = rom.system
print()
print(
f"Modal reduction: {rom.n_full_states} -> {reduced.n_states} states, "
f"energy retained {rom.energy_retained:.4f}"
)
for f, z in zip(rom.frequencies_hz, rom.damping_ratios, strict=True):
print(f" {f:.4f} Hz, damping {z:.4f}")
bundle = write_rom(rom, data / "rom_8ms")
print(f" written: {sorted(p.name for p in bundle.iterdir())}")
# 4. A Kalman filter on the two tower-top accelerometers.
rows = [
i
for i, name in enumerate(reduced.output_names)
if any(tag in name for tag in SENSORS)
]
assert reduced.b is not None and reduced.c is not None and reduced.d is not None
plant = StateSpace(
a=reduced.a,
b=reduced.b,
c=reduced.c[rows],
d=reduced.d[rows],
state_names=reduced.state_names,
input_names=reduced.input_names,
output_names=[reduced.output_names[i] for i in rows],
)
q_c = modal_process_noise([1.0e-4] * len(rom.frequencies_hz))
q_d = discretize_process_noise(plant.a, q_c, DT)
r = measurement_noise_from_specs(SENSOR_STD)
model = build_kalman_model(
plant, CovarianceSet(q=q_d, r=r, p0=np.eye(plant.n_states)), dt=DT
)
assert model.c is not None
# A simulated record from the same model: process and sensor noise only.
rng = np.random.default_rng(1)
truth = np.zeros((N_STEPS, plant.n_states))
state = np.zeros(plant.n_states)
noise_q = np.linalg.cholesky(q_d + 1e-15 * np.eye(plant.n_states))
for k in range(N_STEPS):
state = model.a @ state + noise_q @ rng.standard_normal(plant.n_states)
truth[k] = state
measurements = truth @ model.c.T + rng.standard_normal((N_STEPS, len(rows))) * (
np.asarray(SENSOR_STD)
)
controls = np.zeros((N_STEPS, plant.n_inputs)) # no input perturbation
run = run_filter(model, measurements, controls=controls)
consistency = filter_consistency(run)
tower = reduced.c[[reduced.output_names.index("ED TwrBsMyt, (kN-m)")]]
load, load_std = reconstruct_load(run, tower)
true_load = truth @ tower.T
error = load[:, 0] - true_load[:, 0]
settled = slice(N_STEPS // 10, None)
rmse = float(np.sqrt(np.mean(error[settled] ** 2)))
inside = float(np.mean(np.abs(error[settled]) <= 2.0 * load_std[settled, 0]))
print()
print(f"Kalman filter on {plant.output_names}, dt = {DT} s, {N_STEPS} steps")
print(
f" mean NIS {consistency.mean_nis:.3f} (expected "
f"{consistency.expected_nis:.1f}), consistent: {consistency.consistent}"
)
print(
f" TwrBsMyt: true std {np.std(true_load[settled]):.1f} kN-m, "
f"RMSE {rmse:.1f} kN-m, within 2 sigma {100 * inside:.1f} %"
)
# 5. A mode-tracked family over the sweep, scheduled on wind speed.
systems = [result.state_space(i) for i in range(len(winds))]
family = build_rom_family(
systems,
winds,
n_modes=4,
frequency_band_hz=BAND_HZ,
rank_output="TwrBsMyt",
)
print()
print("Mode-tracked family: frequency (Hz) of each mode")
print(" wind " + "".join(f" mode {k}" for k in range(4)))
for wind, freqs in zip(
family.operating_points, family.frequencies_hz, strict=True
):
print(f" {wind:5.1f} " + "".join(f" {f:6.4f}" for f in freqs))
lowest = np.round(family.tracking_mac.min(axis=0), 4)
print(f" lowest tracking MAC per mode: {lowest}")
scheduled = schedule_model(family, 9.5)
print(
f"Scheduled at 9.5 m/s: {scheduled.system.n_states} states, "
f"min tracking MAC {scheduled.min_tracking_mac:.4f}, "
f"extrapolated {scheduled.extrapolated}"
)
save_rom_family(family, data / "family.npz")
reloaded = load_rom_family(data / "family.npz")
same = np.array_equal(reloaded.systems[3].a, family.systems[3].a)
print(f"Saved and reloaded family.npz: matrices identical {same}")
# 6. A named physical-state family straight from the .lin files.
print()
try:
build_physical_lpv_family(
data / "sweep",
states=["1st tower fore-aft", "Variable speed generator"],
inputs=["horizontal wind speed", "collective blade-pitch"],
outputs=["YawBrTAxp", "GenSpeed"],
)
except ValueError as error:
print("Tower and generator only raises ValueError:")
text = textwrap.shorten(str(error), width=300)
print(
textwrap.fill(
text, width=86, initial_indent=" ", subsequent_indent=" "
)
)
plant_family = build_physical_lpv_family(
data / "sweep",
states=[
"1st tower fore-aft",
"Variable speed generator",
"1st flapwise bending-mode DOF of blade 1",
"1st flapwise bending-mode DOF of blade 2",
"1st flapwise bending-mode DOF of blade 3",
],
inputs=["horizontal wind speed", "collective blade-pitch"],
outputs=["YawBrTAxp", "GenSpeed"],
)
print(
f"Physical family: {len(plant_family.systems)} points, "
f"{len(plant_family.state_names)} states, scheduled on "
f"{plant_family.scheduling_variable} ({plant_family.scheduling_unit})"
)
print(f" twins added: {len(plant_family.auto_added_states)}")
print(
" dropped coupling: "
f"max {float(np.max(plant_family.dropped_coupling)):.2e}"
)
save_lpv_family(plant_family, data / "physical_family.npz")
reloaded_plant = load_lpv_family(data / "physical_family.npz")
same = reloaded_plant.state_names == plant_family.state_names and all(
np.array_equal(getattr(x, m), getattr(y, m))
for x, y in zip(reloaded_plant.systems, plant_family.systems, strict=True)
for m in ("a", "b", "c", "d")
)
print(f" saved and reloaded physical_family.npz: matrices identical {same}")
with figure_style():
figure, axes = plt.subplots(figsize=(6.5, 4.0))
time = DT * np.arange(N_STEPS)
window = (time >= 200.0) & (time <= 260.0)
axes.plot(time[window], true_load[window, 0], label="true")
axes.plot(time[window], load[window, 0], label="estimate")
axes.fill_between(
time[window],
load[window, 0] - 2.0 * load_std[window, 0],
load[window, 0] + 2.0 * load_std[window, 0],
alpha=0.25,
linewidth=0.0,
label="estimate +/- 2 sigma",
)
axes.set_xlabel("Time (s)")
axes.set_ylabel("TwrBsMyt perturbation (kN-m)")
axes.legend(loc="upper right", fontsize=8)
figure.tight_layout()
save_figure(figure, figures / "tutorial_kalman_load.png")
return 0
Run it
$ python examples/06_state_space_export.py --output export
State space at 8.0 m/s: 12 states, 2 inputs, 8 outputs, stable: False
inputs: ['IfW Extended input: horizontal wind speed (steady/uniform wind) (hub), m/s', 'ED Extended input: collective blade-pitch command, rad']
outputs: ['ED YawBrTAxp, (m/s^2)', 'ED YawBrTAyp, (m/s^2)', 'ED GenSpeed, (rpm)', 'ED TwrBsMyt, (kN-m)', 'ED TwrBsMxt, (kN-m)', 'ED RootMyc1, (kN-m)', 'ED RootMyc2, (kN-m)', 'ED RootMyc3, (kN-m)']
Modal reduction: 12 -> 8 states, energy retained 1.0000
0.3215 Hz, damping 0.0589
0.7133 Hz, damping 0.1089
0.3000 Hz, damping 0.0097
0.5564 Hz, damping 0.1415
written: ['manifest.json', 'modes.csv', 'state_space.mat']
Kalman filter on ['ED YawBrTAxp, (m/s^2)', 'ED YawBrTAyp, (m/s^2)'], dt = 0.05 s, 6000 steps
mean NIS 1.998 (expected 2.0), consistent: True
TwrBsMyt: true std 4422.2 kN-m, RMSE 799.9 kN-m, within 2 sigma 95.9 %
Mode-tracked family: frequency (Hz) of each mode
wind mode 0 mode 1 mode 2 mode 3
4.0 0.3213 0.6936 0.3001 0.6117
5.0 0.3213 0.6973 0.3001 0.5966
6.0 0.3214 0.7019 0.3001 0.5824
7.0 0.3215 0.7072 0.3000 0.5690
8.0 0.3215 0.7133 0.3000 0.5564
9.0 0.3216 0.7202 0.2999 0.5447
10.0 0.3217 0.7278 0.2999 0.5337
11.0 0.3217 0.7324 0.2999 0.5278
13.0 0.3217 0.7324 0.2999 0.5280
15.0 0.3218 0.7324 0.2999 0.5283
lowest tracking MAC per mode: [1. 1. 1. 0.9999]
Scheduled at 9.5 m/s: 8 states, min tracking MAC 1.0000, extrapolated False
Saved and reloaded family.npz: matrices identical True
Tower and generator only raises ValueError:
the state selection drops 0.471 relative coupling at WindSpeed = 11 m/s (A: 0.471,
C: 0.0027), above max_dropped_coupling=0.05; the largest contributing dropped states
are 'ED 1st flapwise bending-mode DOF of blade 1 (internal DOF index = DOF_BF(1,1))'
(0.272), 'ED 1st flapwise bending-mode [...]
Physical family: 10 points, 10 states, scheduled on WindSpeed (m/s)
twins added: 5
dropped coupling: max 1.56e-18
saved and reloaded physical_family.npz: matrices identical True
Figure written: tutorial_kalman_load.png
The tower-base fore-aft moment reconstructed from two tower-top accelerometers: the simulated truth, the filter estimate, and the estimate’s two-standard-deviation band, over 60 s of the 300 s record.
Step by step
1. The sweep. The script writes the ten operating points and runs
ModalPipeline on the wind-speed axis, because the rotor speed is constant
above rated (see Tutorial 3: Mode identification and the Campbell diagram).
2. The full-order model. state_space() returns
the azimuth-averaged model of one operating point as a
StateSpace, with the .lin names of its states,
inputs and outputs. The three blade-root moments are in multi-blade coordinates
(collective, cosine, sine), as recorded in system.output_mbc_coordinates. The model
is not asymptotically stable (stable: False) because the generator-azimuth state
has a zero eigenvalue: the rotor position is free. A filter should not run on such a
model, which is the reason for the next step.
3. Modal reduction. modal_reduction() keeps the underdamped
modes in the band 0.1 to 2 Hz with the largest residue in the tower-base fore-aft
moment TwrBsMyt, as a real block-diagonal model with two states per mode and the
original inputs and outputs. Four modes are kept: the two tower modes and the
regressive and collective flap modes; the progressive mode, which drives the tower
base least, is left out. The azimuth and rotor-speed states are gone.
write_rom() writes the matrices (state_space.mat), a table of
the retained modes (modes.csv) and a manifest with the VANE version.
4. A Kalman filter. The sensors are the two tower-top accelerometers,
YawBrTAxp and YawBrTAyp. The script keeps their rows of the reduced model,
builds the continuous process noise of each mode with
modal_process_noise(), discretizes it exactly over the
0.05 s step with discretize_process_noise(), takes the
sensor noise (0.02 m/s²) with
measurement_noise_from_specs(), and assembles the
discrete filter model with build_kalman_model().
A 300 s record is then simulated from the same model, with process and sensor noise.
run_filter() runs the filter over it (the inputs are
zero perturbations, passed explicitly because the accelerometers have a feedthrough
term). filter_consistency() tests the normalized
innovation squared: its mean of 1.998 against the expected 2.0 (two sensors) shows the
filter is tuned to the noise it sees.
reconstruct_load() reconstructs TwrBsMyt, which no
sensor measures, from the filter state, with its standard deviation. The error is 800
kN-m against a load of 4422 kN-m standard deviation, and 95.9 % of the errors lie
within two standard deviations, close to the 95.4 % expected of a Gaussian error.
Because the truth comes from the filter’s own model, this is a check of the
implementation, not of the accuracy to expect on a real turbine.
5. A mode-tracked family. build_rom_family() reduces
every operating point to the same four modes, matched by their output residues, in the
same state order and phase, so that neighbouring members can be blended. The tracking
MAC of every mode is 1.0 (0.9999 for the regressive mode). The flap frequencies change
with rotor speed up to 11 m/s and stay constant above rated. Mode 3 is the regressive
mode, whose frequency falls from 0.612 to 0.528 Hz.
schedule_model() interpolates the family at 9.5 m/s between its
9 and 10 m/s members; min_tracking_mac and extrapolated tell whether the blend
can be trusted. save_rom_family() stores the family in a
.npz file that load_rom_family() reads back with
allow_pickle=False.
6. A physical-state family. A controller usually needs named states rather than
modal ones. vane.export.build_physical_lpv_family() reads the .lin directory,
averages each operating point’s azimuth files without the MBC3 transform, and selects
the named states, inputs and outputs by substring. Selecting only the tower fore-aft
and generator states is refused: the flap states that are left out carry 47 % of the
coupling into the selected rows, far above the default ceiling of 5 %
(max_dropped_coupling=0.05), and the error names the states responsible. Adding the
three flap states makes the selection complete: the side-side tower state, still left
out, couples to the blades through the sine of their azimuth, which averages to zero
over the twelve azimuths, so the dropped coupling falls to rounding (1.6e-18).
pair_derivatives adds the rate twin of each requested displacement, so the five
requested states become ten. save_lpv_family() stores the
LPVFamily in a .npz file, with its channel names,
scheduling variable and unit, and the schema version of the archive;
load_lpv_family() reads it back with allow_pickle=False.
Exercises
Rank the reduction by the blade-root moment instead (
rank_output="RootMyc1", the collective coordinate after the transform). The collective flap mode now comes first, and the progressive mode replaces the regressive one.Lower
SENSOR_STDto 0.005 m/s². The filter stays consistent and the error of the reconstructed moment falls from 800 to 635 kN-m: the process noise, not the sensors, now limits the estimate.
See State-space export for the options of every step, and Realization, reduction, and observability and LPV export, load maps, and estimation for the methods.