Linear-quadratic regulated oscillator
A harmonic oscillator that a linear-quadratic regulator brings to rest.
- Name:
lqr_oscillator - Field:
engineering/control - Formalism: ordinary differential equations
- References: Åström and Murray 2008
Factory
@register_model(
field='engineering/control',
tags=('control', 'feedback', 'linear-quadratic regulator', 'oscillator', 'ode'),
references=('Åström and Murray 2008',),
defaults={'until': 10.0, 'dt': 0.05, 'solver': RK4Solver(dt=0.001)},
observables=lambda system: [ObservableSpec('energy', _oscillator_energy)],
)
def lqr_oscillator(
x: float = 1.0,
*,
omega: float = 1.0,
zeta: float = 0.0,
q: float = 1.0,
r: float = 0.1,
) -> System:
"""Return a harmonic oscillator that a linear-quadratic regulator brings to rest.
The gain `K` of the regulator minimises `∫ (q (x^2 + x'^2) + r force^2) dt` for the linearisation of the
oscillator, which `LinearizationAnalysis` computes from the Riccati equation, and the state feedback sets
`force = -K (x, x')`. A smaller weight `r` of the force gives a larger gain and a faster return to rest.
Args:
x:
The initial displacement.
omega:
The natural frequency.
zeta:
The damping ratio without control.
q:
The weight of the displacement and the velocity in the cost.
r:
The weight of the force in the cost.
Returns:
system:
The model, with the displacement `oscillator.x`, the force `oscillator.force`, and the observable
`energy`.
Examples:
>>> from skcomplex.models.engineering.control import lqr_oscillator
>>> [round(value, 3) for value in lqr_oscillator().interactions[0].gain]
[2.317, 3.825]
"""
oscillator = DrivenOscillatorElement('oscillator', x=x, omega=omega, zeta=zeta)
model = LinearizationAnalysis(['oscillator.force'], state_weights=q * np.eye(2), control_weights=[[r]])
model.check().prepare(System(elements=[oscillator])).run()
gain = model.gain_
if gain is None:
error_msg = 'The linearization of the oscillator gives no gain, since it has no cost weights.'
raise ValidationError(error_msg)
states = tuple(name.removeprefix('oscillator.') for name in model.states_)
regulator = StateFeedbackInteraction(
'regulator',
plant='oscillator',
control='force',
states=states,
gain=tuple(float(value) for value in gain[0]),
limits=(-FORCE_LIMIT, FORCE_LIMIT),
)
return System(elements=[oscillator], interactions=[regulator])
Components
class DrivenOscillatorElement(Element):
"""Damped harmonic oscillator that a bounded force drives, `x'' = -omega^2 x - 2 zeta omega x' + force`.
Args:
label:
The label of the oscillator.
visual:
How the element is drawn, or `None` for `default_visual_`.
x:
The initial displacement.
x_d1:
The initial velocity.
force:
The force per unit mass, the control input.
omega:
The natural frequency.
zeta:
The damping ratio.
Examples:
>>> from skcomplex.models.engineering.control import DrivenOscillatorElement
>>> DrivenOscillatorElement('oscillator', x=1.0).omega
1.0
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
x: Annotated[float, StateSpec(order=2)] = 0.0,
x_d1: float = 0.0,
force: Annotated[float, ParameterSpec(adaptive=True, bounds=(-FORCE_LIMIT, FORCE_LIMIT))] = 0.0,
omega: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
zeta: Annotated[float, ParameterSpec(domain=RealSet.closed_open(0, oo))] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.x = x
self.x_d1 = x_d1
self.force = force
self.omega = omega
self.zeta = zeta
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the rate of the velocity.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
return [RateLaw(s.x, -(s.omega**2) * s.x - 2 * s.zeta * s.omega * derivative(s.x, 1) + s.force)]
class StateFeedbackInteraction(Interaction):
"""Linear state feedback `u = -K x` that sets a control input of a plant from its states.
The control input is an adaptive parameter of the plant, which the controller sets after every solver step, so it
holds within a step like a sampled controller. A gain from `LinearizationAnalysis` of the plant makes the
controller a linear-quadratic regulator.
Args:
label:
The label of the controller.
plant:
The path of the element that holds the states and the control input.
control:
The name of the adaptive parameter of the plant that the controller sets.
states:
The names of the fed back states of the plant, such as `('x', 'x_d1')` for a position and its velocity.
gain:
The gain of every state, in the order of `states`.
limits:
The lower and upper limit of the control input, or `None` for no limits.
Examples:
>>> from skcomplex.models.engineering.control import StateFeedbackInteraction
>>> StateFeedbackInteraction(
... 'regulator', plant='mass', control='force', states=('x', 'x_d1'), gain=(1.0, 1.7)
... ).gain
(1.0, 1.7)
"""
def __init__(
self,
label: str,
*,
plant: Annotated[str, RoleSpec()],
control: str = 'u',
states: Sequence[str] = (),
gain: Sequence[float] = (),
limits: tuple[float, float] | None = None,
) -> None:
super().__init__(label)
self.plant = plant
self.control = control
self.states = states
self.gain = gain
self.limits = limits
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the control input the feedback sets.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
feedback = -sum(
(float(gain) * getattr(s.plant, name) for name, gain in zip(self.states, self.gain, strict=True)),
0,
)
return [BoundaryLaw(getattr(s.plant, self.control), clip(feedback, *(self.limits or (None, None))))]
Equations
oscillator.x: diff(x_oscillator(t), t, t) == -omega_oscillator^2*x_oscillator(t) - 2*omega_oscillator*zeta_oscillator*diff(x_oscillator(t), t) + force_oscillator(t)
oscillator.force: diff(force_oscillator(t), t) == 0