Control

Control: feedback loops that keep or drive a system towards a goal, and the controllers that close them.

BangBangPlantElement(label, visual=None, *, x=0.0, y=0.0, u=1.0, a=0.0)

Bases: Element

Second-order linear plant under a bang-bang control, held between samples, that steers it to the origin.

The plant follows x' = y and y' = a x + u with the control u of size 1, a = 0 for a free particle, a = -1 for a harmonic oscillator, and a = 1 for a saddle. The control stays constant between samples, and at every sample takes -1 above the switching curve y = g(x) of the plant and 1 below it. The curve is made of the two trajectories of the plant under u = 1 and u = -1 that reach the origin: g(x) = -sign(x) sqrt(2 |x|) for the free particle, the half circles g(x) = sqrt(1 - (x + 1)^2) on [-2, 0] and -sqrt(1 - (x - 1)^2) on [0, 2], and 0 outside, for the harmonic oscillator, and g(x) = sqrt(x^2 - 2 x) for x <= 0 and -sqrt(x^2 + 2 x) for x >= 0 for the saddle.

Parameters:

Name	Type	Description	Default
label	str	The label of the plant.	required
visual	VisualSpec | None	How the plant is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The initial position.	0.0
y	Annotated[float, StateSpec()]	The initial velocity.	0.0
u	Annotated[float, StateSpec(domain=FiniteEnumeratedSet([-1, 1]))]	The initial control, -1 or 1.	1.0
a	Annotated[float, ParameterSpec(domain=FiniteEnumeratedSet([-1, 0, 1]))]	The coefficient of the position in the acceleration, 0, -1, or 1.	0.0

Examples:

>>> from skcomplex.systems.engineering.control import BangBangPlantElement
>>> BangBangPlantElement('plant', a=-1.0).a
-1.0

Source code in src/skcomplex/systems/engineering/control/_minimum_time_control.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    y: Annotated[float, StateSpec()] = 0.0,
    u: Annotated[float, StateSpec(domain=FiniteEnumeratedSet([-1, 1]))] = 1.0,
    a: Annotated[float, ParameterSpec(domain=FiniteEnumeratedSet([-1, 0, 1]))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y
    self.u = u
    self.a = a

CartPendulumElement(label, visual=None, *, x=0.0, x_d1=0.0, theta=0.0, theta_d1=0.0, force=0.0, M=1.0, m=0.1, l=0.5, g=GRAVITY)

Bases: Element

Inverted pendulum on a cart that a bounded horizontal force pushes, after Åström and Murray.

The cart of mass M at position x carries a point mass m at the end of a rod of length l, at the angle theta from upright, positive towards negative x. The equations of motion are (M + m) x'' - m l theta'' cos(theta) + m l theta'^2 sin(theta) = force and l theta'' - x'' cos(theta) = g sin(theta), so upright is unstable and a positive force tilts the pendulum back from a positive angle.

Parameters:

Name	Type	Description	Default
label	str	The label of the cart.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(order=2)]	The initial position of the cart.	0.0
x_d1	float	The initial velocity of the cart.	0.0
theta	Annotated[float, StateSpec(order=2)]	The initial angle from upright.	0.0
theta_d1	float	The initial angular velocity.	0.0
force	Annotated[float, ParameterSpec(adaptive=True, bounds=(-FORCE_LIMIT, FORCE_LIMIT))]	The horizontal force on the cart, the control input.	0.0
M	Annotated[float, ParameterSpec(domain=open(0, oo))]	The mass of the cart.	1.0
m	Annotated[float, ParameterSpec(domain=open(0, oo))]	The mass at the end of the rod.	0.1
l	Annotated[float, ParameterSpec(domain=open(0, oo))]	The length of the rod.	0.5
g	Annotated[float, ParameterSpec(domain=open(0, oo))]	The acceleration of gravity.	GRAVITY

Examples:

>>> from skcomplex.systems.engineering.control import CartPendulumElement
>>> CartPendulumElement('cart', theta=0.1).theta
0.1

Source code in src/skcomplex/systems/engineering/control/_regulation.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec(order=2)] = 0.0,
    x_d1: float = 0.0,
    theta: Annotated[float, StateSpec(order=2)] = 0.0,
    theta_d1: float = 0.0,
    force: Annotated[float, ParameterSpec(adaptive=True, bounds=(-FORCE_LIMIT, FORCE_LIMIT))] = 0.0,
    M: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 1.0,
    m: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.1,
    l: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = 0.5,
    g: Annotated[float, ParameterSpec(domain=RealSet.open(0, oo))] = GRAVITY,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.x_d1 = x_d1
    self.theta = theta
    self.theta_d1 = theta_d1
    self.force = force
    self.M = M
    self.m = m
    self.l = l
    self.g = g

ControlInputElement(label, visual=None, *, u=0.0)

Bases: Element

Control input that a plant reads, a bounded adaptive parameter scaled to the interval from -1 to 1.

Parameters:

Name	Type	Description	Default
label	str	The label of the input.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
u	Annotated[float, ParameterSpec(adaptive=True, bounds=(-1.0, 1.0))]	The value of the input.	0.0

Examples:

>>> from skcomplex.systems.engineering.control import ControlInputElement
>>> ControlInputElement('input', u=0.5).u
0.5

Source code in src/skcomplex/systems/engineering/control/_koopman.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    u: Annotated[float, ParameterSpec(adaptive=True, bounds=(-1.0, 1.0))] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.u = u

ConvectionElement(label, visual=None, *, a_T=0.0, a_Q=0.0, h=1.0, fluid=0.0)

Bases: Element

Surface that exchanges heat with a fluid of a fixed temperature, with a_Q = h (a_T - fluid).

Parameters:

Name	Type	Description	Default
label	str	The label of the element, unique among the parts of its system.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
a_T	Annotated[object, StateSpec(algebraic=True)]	The temperature of the port.	0.0
a_Q	Annotated[object, StateSpec(algebraic=True)]	The heat flow into the port.	0.0
h	Annotated[object, ParameterSpec(domain=POSITIVE)]	The convective conductance of the surface.	1.0
fluid	Annotated[object, ParameterSpec()]	The temperature of the fluid.	0.0

Examples:

>>> from skcomplex.systems.engineering.control import ConvectionElement
>>> ConvectionElement('outside', fluid=5.0).fluid
5.0

Source code in src/skcomplex/systems/engineering/control/_control.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    a_T: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    a_Q: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    h: Annotated[object, ParameterSpec(domain=POSITIVE)] = 1.0,
    fluid: Annotated[object, ParameterSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.a_T = a_T
    self.a_Q = a_Q
    self.h = h
    self.fluid = fluid

DrivenOscillatorElement(label, visual=None, *, x=0.0, x_d1=0.0, force=0.0, omega=1.0, zeta=0.0)

Bases: Element

Damped harmonic oscillator that a bounded force drives, x'' = -omega^2 x - 2 zeta omega x' + force.

Parameters:

Name	Type	Description	Default
label	str	The label of the oscillator.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec(order=2)]	The initial displacement.	0.0
x_d1	float	The initial velocity.	0.0
force	Annotated[float, ParameterSpec(adaptive=True, bounds=(-FORCE_LIMIT, FORCE_LIMIT))]	The force per unit mass, the control input.	0.0
omega	Annotated[float, ParameterSpec(domain=open(0, oo))]	The natural frequency.	1.0
zeta	Annotated[float, ParameterSpec(domain=closed_open(0, oo))]	The damping ratio.	0.0

Examples:

>>> from skcomplex.systems.engineering.control import DrivenOscillatorElement
>>> DrivenOscillatorElement('oscillator', x=1.0).omega
1.0

Source code in src/skcomplex/systems/engineering/control/_regulation.py

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

HeatCapacityElement(label, visual=None, *, a_T=0.0, a_Q=0.0, T=20.0, C=1.0)

Bases: Element

Body of one temperature T, with C T' = a_Q and T = a_T.

Parameters:

Name	Type	Description	Default
label	str	The label of the element, unique among the parts of its system.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
a_T	Annotated[object, StateSpec(algebraic=True)]	The temperature of the port.	0.0
a_Q	Annotated[object, StateSpec(algebraic=True)]	The heat flow into the port.	0.0
T	Annotated[object, StateSpec()]	The initial temperature of the body.	20.0
C	Annotated[object, ParameterSpec(domain=POSITIVE)]	The heat capacity.	1.0

Examples:

>>> from skcomplex.systems.engineering.control import HeatCapacityElement
>>> HeatCapacityElement('room', T=15.0).T
15.0

Source code in src/skcomplex/systems/engineering/control/_control.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    a_T: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    a_Q: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    T: Annotated[object, StateSpec()] = 20.0,
    C: Annotated[object, ParameterSpec(domain=POSITIVE)] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.a_T = a_T
    self.a_Q = a_Q
    self.T = T
    self.C = C

HeatingResistorElement(label, visual=None, *, p_v=0.0, p_i=0.0, n_v=0.0, n_i=0.0, h_T=0.0, h_Q=0.0, on=1.0, R=1.0)

Bases: TwoPinElement

Resistor that delivers its electric power to a thermal port h while its switch on is closed.

The current is on (p_v - n_v) / R, and the heat flow into the thermal port is -(p_v - n_v) p_i.

Parameters:

Name	Type	Description	Default
label	str	The label of the element, unique among the parts of its system.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
p_v	object	The potential of the positive pin.	0.0
p_i	object	The current into the positive pin.	0.0
n_v	object	The potential of the negative pin.	0.0
n_i	object	The current into the negative pin.	0.0
h_T	Annotated[object, StateSpec(algebraic=True)]	The temperature of the thermal port.	0.0
h_Q	Annotated[object, StateSpec(algebraic=True)]	The heat flow into the thermal port.	0.0
on	Annotated[object, StateSpec()]	The initial position of the switch, 1 for closed and 0 for open.	1.0
R	Annotated[object, ParameterSpec(domain=POSITIVE)]	The resistance.	1.0

Examples:

>>> from skcomplex.systems.engineering.control import HeatingResistorElement
>>> sorted(HeatingResistorElement.ports)
['h', 'n', 'p']

Source code in src/skcomplex/systems/engineering/control/_control.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    p_v: object = 0.0,
    p_i: object = 0.0,
    n_v: object = 0.0,
    n_i: object = 0.0,
    h_T: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    h_Q: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    on: Annotated[object, StateSpec()] = 1.0,
    R: Annotated[object, ParameterSpec(domain=POSITIVE)] = 1.0,
) -> None:
    super().__init__(label, visual, p_v=p_v, p_i=p_i, n_v=n_v, n_i=n_i)
    self.h_T = h_T
    self.h_Q = h_Q
    self.on = on
    self.R = R

IntegratorElement(label, visual=None, *, value=0.0)

Bases: Element

Element that holds the integral of a signal, whose rate a controller sets.

Parameters:

Name	Type	Description	Default
label	str	The label of the integrator.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
value	Annotated[float, StateSpec()]	The initial value of the integral.	0.0

Examples:

>>> from skcomplex.systems.engineering.control import IntegratorElement
>>> IntegratorElement('memory').value
0.0

Source code in src/skcomplex/systems/engineering/control/_regulation.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    value: Annotated[float, StateSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.value = value

PIDInteraction(label, *, plant, integrator, kp=1.0, ki=0.0, kd=0.0, setpoint=0.0, control='u', measure='y', limits=None)

Bases: Interaction

Proportional-integral-derivative controller that sets a control input of a plant from the error of a measure.

The error is e = setpoint - y for the measured state y, and the controller sets the control input to u = kp e + ki I - kd y', clipped to the limits, with the derivative taken on the measure. The integral I' = e of the error lives in an IntegratorElement bound to the role integrator. 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. The derivative term reads the first derivative of the measure, so a measure with a derivative term is a state of order 2.

Parameters:

Name	Type	Description	Default
label	str	The label of the controller.	required
plant	Annotated[str, RoleSpec()]	The path of the element that holds the measure and the control input.	required
integrator	Annotated[str, RoleSpec()]	The path of the integrator of the error.	required
kp	Annotated[float, ParameterSpec()]	The proportional gain.	1.0
ki	Annotated[float, ParameterSpec()]	The integral gain.	0.0
kd	Annotated[float, ParameterSpec()]	The derivative gain.	0.0
setpoint	Annotated[float, ParameterSpec()]	The value the measure is kept at.	0.0
control	str	The name of the adaptive parameter of the plant that the controller sets.	'u'
measure	str	The name of the measured state of the plant.	'y'
limits	tuple[float, float] | None	The lower and upper limit of the control input, or None for no limits.	None

Examples:

>>> from skcomplex.systems.engineering.control import PIDInteraction
>>> PIDInteraction('pid', plant='cart', integrator='memory', control='force', measure='theta', kp=40.0).kp
40.0

Source code in src/skcomplex/systems/engineering/control/_regulation.py

def __init__(
    self,
    label: str,
    *,
    plant: Annotated[str, RoleSpec()],
    integrator: Annotated[str, RoleSpec()],
    kp: Annotated[float, ParameterSpec()] = 1.0,
    ki: Annotated[float, ParameterSpec()] = 0.0,
    kd: Annotated[float, ParameterSpec()] = 0.0,
    setpoint: Annotated[float, ParameterSpec()] = 0.0,
    control: str = 'u',
    measure: str = 'y',
    limits: tuple[float, float] | None = None,
) -> None:
    super().__init__(label)
    self.plant = plant
    self.integrator = integrator
    self.kp = kp
    self.ki = ki
    self.kd = kd
    self.setpoint = setpoint
    self.control = control
    self.measure = measure
    self.limits = limits

PlantStateElement(label, visual=None, *, x=0.0, y=0.0)

Bases: Element

State of a plant with two variables, whose rates a pair of interactions gives.

Parameters:

Name	Type	Description	Default
label	str	The label of the plant.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
x	Annotated[float, StateSpec()]	The value of the first variable.	0.0
y	Annotated[float, StateSpec()]	The value of the second variable, the output that tracks the reference.	0.0

Examples:

>>> from skcomplex.systems.engineering.control import PlantStateElement
>>> PlantStateElement('plant', x=0.5).x
0.5

Source code in src/skcomplex/systems/engineering/control/_koopman.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    x: Annotated[float, StateSpec()] = 0.0,
    y: Annotated[float, StateSpec()] = 0.0,
) -> None:
    super().__init__(label, visual=visual)
    self.x = x
    self.y = y

StateFeedbackInteraction(label, *, plant, control='u', states=(), gain=(), limits=None)

Bases: 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.

Parameters:

Name	Type	Description	Default
label	str	The label of the controller.	required
plant	Annotated[str, RoleSpec()]	The path of the element that holds the states and the control input.	required
control	str	The name of the adaptive parameter of the plant that the controller sets.	'u'
states	Sequence[str]	The names of the fed back states of the plant, such as ('x', 'x_d1') for a position and its velocity.	()
gain	Sequence[float]	The gain of every state, in the order of states.	()
limits	tuple[float, float] | None	The lower and upper limit of the control input, or None for no limits.	None

Examples:

>>> from skcomplex.systems.engineering.control import StateFeedbackInteraction
>>> StateFeedbackInteraction(
...     'regulator', plant='mass', control='force', states=('x', 'x_d1'), gain=(1.0, 1.7)
... ).gain
(1.0, 1.7)

Source code in src/skcomplex/systems/engineering/control/_regulation.py

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

ThermalConductorElement(label, visual=None, *, a_T=0.0, a_Q=0.0, b_T=0.0, b_Q=0.0, G=1.0)

Bases: Element

Conductor between two ports, with the heat flow a_Q = G (a_T - b_T) and a_Q + b_Q = 0.

Parameters:

Name	Type	Description	Default
label	str	The label of the element, unique among the parts of its system.	required
visual	VisualSpec | None	How the element is drawn, or None for default_visual_.	None
a_T	Annotated[object, StateSpec(algebraic=True)]	The temperature of the first port.	0.0
a_Q	Annotated[object, StateSpec(algebraic=True)]	The heat flow into the first port.	0.0
b_T	Annotated[object, StateSpec(algebraic=True)]	The temperature of the second port.	0.0
b_Q	Annotated[object, StateSpec(algebraic=True)]	The heat flow into the second port.	0.0
G	Annotated[object, ParameterSpec(domain=POSITIVE)]	The thermal conductance.	1.0

Examples:

>>> from skcomplex.systems.engineering.control import ThermalConductorElement
>>> ThermalConductorElement('wall', G=2.0).G
2.0

Source code in src/skcomplex/systems/engineering/control/_control.py

def __init__(
    self,
    label: str,
    visual: VisualSpec | None = None,
    *,
    a_T: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    a_Q: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    b_T: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    b_Q: Annotated[object, StateSpec(algebraic=True)] = 0.0,
    G: Annotated[object, ParameterSpec(domain=POSITIVE)] = 1.0,
) -> None:
    super().__init__(label, visual=visual)
    self.a_T = a_T
    self.a_Q = a_Q
    self.b_T = b_T
    self.b_Q = b_Q
    self.G = G

heated_room(T=10.0, *, setting=20.0, band=1.0, outside=5.0, V=230.0, R=25.0, C=200000.0, wall_capacity=1000000.0, wall_conductance=100.0, surface_conductance=50.0)

Return a room assembled from thermal and electrical components, heated under a thermostat with hysteresis.

A voltage source drives a heating resistor whose heat enters the room, and the room loses heat through a wall of its own heat capacity to the outside air. The thermostat opens the switch of the resistor when the room rises through setting + band / 2 and closes it when the room falls through setting - band / 2.

Parameters:

Name	Type	Description	Default
T	float	The initial temperature of the room.	10.0
setting	float	The temperature the thermostat keeps.	20.0
band	float	The width of the hysteresis band around the setting.	1.0
outside	float	The temperature of the outside air, which is also the initial temperature of the wall.	5.0
V	float	The voltage of the source.	230.0
R	float	The resistance of the heater.	25.0
C	float	The heat capacity of the room.	200000.0
wall_capacity	float	The heat capacity of the wall.	1000000.0
wall_conductance	float	The thermal conductance between the room and the wall.	100.0
surface_conductance	float	The convective conductance between the wall and the outside air.	50.0

Returns:

Name	Type	Description
system	System	The system, with the room temperature room.T and the switch heater.on, in seconds.

Examples:

>>> from skcomplex.systems.engineering.control import heated_room
>>> heated_room(setting=22.0).interactions[-1].direction
-1

Source code in src/skcomplex/systems/engineering/control/_control.py

@register_system(
    field='engineering/control',
    tags=('control', 'thermal', 'acausal components', 'hysteresis', 'ode', 'events'),
    references=('Fritzson 2015',),
    defaults={'until': 14400.0, 'dt': 10.0},
    observables=lambda system: [ObservableSpec('power', _heating_power)],
)
def heated_room(
    T: float = 10.0,
    *,
    setting: float = 20.0,
    band: float = 1.0,
    outside: float = 5.0,
    V: float = 230.0,
    R: float = 25.0,
    C: float = 2e5,
    wall_capacity: float = 1e6,
    wall_conductance: float = 100.0,
    surface_conductance: float = 50.0,
) -> System:
    """Return a room assembled from thermal and electrical components, heated under a thermostat with hysteresis.

    A voltage source drives a heating resistor whose heat enters the room, and the room loses heat through a wall of its
    own heat capacity to the outside air. The thermostat opens the switch of the resistor when the room rises through
    `setting + band / 2` and closes it when the room falls through `setting - band / 2`.

    Args:
        T:
            The initial temperature of the room.

        setting:
            The temperature the thermostat keeps.

        band:
            The width of the hysteresis band around the setting.

        outside:
            The temperature of the outside air, which is also the initial temperature of the wall.

        V:
            The voltage of the source.

        R:
            The resistance of the heater.

        C:
            The heat capacity of the room.

        wall_capacity:
            The heat capacity of the wall.

        wall_conductance:
            The thermal conductance between the room and the wall.

        surface_conductance:
            The convective conductance between the wall and the outside air.

    Returns:
        system:
            The system, with the room temperature `room.T` and the switch `heater.on`, in seconds.

    Examples:
        >>> from skcomplex.systems.engineering.control import heated_room
        >>> heated_room(setting=22.0).interactions[-1].direction
        -1
    """
    upper, lower = setting + band / 2, setting - band / 2
    return System(
        elements=[
            VoltageSourceElement('source', V=V),
            HeatingResistorElement('heater', R=R),
            GroundElement('ground'),
            HeatCapacityElement('room', T=T, C=C),
            ThermalConductorElement('inner', G=wall_conductance),
            HeatCapacityElement('wall', T=outside, C=wall_capacity),
            ConvectionElement('surface', h=surface_conductance, fluid=outside),
        ],
        interactions=[
            ConnectionInteraction('supply', members=['source', 'heater'], ports=['p', 'p']),
            ConnectionInteraction('neutral', members=['source', 'heater', 'ground'], ports=['n', 'n', 'p']),
            ConnectionInteraction('air', members=['heater', 'room', 'inner'], ports=['h', 'a', 'a']),
            ConnectionInteraction('masonry', members=['inner', 'wall', 'surface'], ports=['b', 'a', 'a']),
            CrossingInteraction(
                'open_switch',
                target='heater',
                guard=lambda s: s.target.h_T - upper,
                reset={'on': lambda s: 0.0},
                direction=1,
            ),
            CrossingInteraction(
                'close_switch',
                target='heater',
                guard=lambda s: s.target.h_T - lower,
                reset={'on': lambda s: 1.0},
                direction=-1,
            ),
        ],
    )

koopman_model_predictive_control(plant='duffing', *, controlled=True, dictionary=None, x=0.5, y=0.0, n_runs=100, duration=1000, random_state=0)

Return a plant that a discrete-time controller steers with model predictive control on a fitted surrogate.

The factory runs the plant under a random binary input at the limits of its range, drawn every sample, and fits ModeDecompositionAnalysis with the input to the runs. The controller lifts the states with the dictionary, condenses the quadratic tracking cost over the horizon with the matrices of the surrogate, solves the bounded least squares problem with scipy.optimize.lsq_linear, and sets the input to the first value of the solution. The plant reads the input through a role, and CoSimulationSolver holds it between samples. The plants are those of figure 6.6 of the review of Brunton et al.: the forced Duffing oscillator x'' = a x + b x^3 + d x' + f cos(w t) + (3 + cos(x)) u, the Van der Pol oscillator x'' = -x - mu (x^2 - 1) x' + u, and the bilinear DC motor, whose second variable is the output that tracks the reference. The motor reads four times the input, as Korda and Mezić scale it. One time unit of the system is one sample of the controller.

Parameters:

Name	Type	Description	Default
plant	str	'duffing', 'van_der_pol', or 'motor'.	'duffing'
controlled	bool	Whether the controller solves the tracking problem, or holds the input at 0 with no fit.	True
dictionary	TransformerMixin | None	A scikit-learn transformer that computes the observables and keeps the variables among them, or None for the monomials up to degree 3.	None
x	float	The initial value of the first variable.	0.5
y	float	The initial value of the second variable.	0.0
n_runs	int	The number of training runs.	100
duration	int	The number of samples of every training run.	1000
random_state	int | Generator | None	The seed or generator of the training runs.	0

Returns:

Name	Type	Description
system	System	The system, labelled by the plant, with the variables plant/state.x and plant/state.y, the input controller/input.u, and the observable reference.

Raises:

Type	Description
ValidationError	The plant is unknown.

Examples:

>>> from skcomplex.systems.engineering.control import koopman_model_predictive_control
>>> system = koopman_model_predictive_control('motor', controlled=False)
>>> [subsystem.label for subsystem in system.systems]
['controller', 'plant']

Source code in src/skcomplex/systems/engineering/control/_koopman.py

@register_system(
    field='engineering/control',
    tags=('control', 'model predictive control', 'koopman', 'co-simulation', 'surrogate'),
    references=('Korda and Mezić 2018', 'Brunton et al. 2022'),
    defaults={'until': 100.0, 'dt': 1.0},
    observables=_observables,
)
def koopman_model_predictive_control(
    plant: str = 'duffing',
    *,
    controlled: bool = True,
    dictionary: TransformerMixin | None = None,
    x: float = 0.5,
    y: float = 0.0,
    n_runs: int = 100,
    duration: int = 1000,
    random_state: int | np.random.Generator | None = 0,
) -> System:
    """Return a plant that a discrete-time controller steers with model predictive control on a fitted surrogate.

    The factory runs the plant under a random binary input at the limits of its range, drawn every sample, and fits
    `ModeDecompositionAnalysis` with the input to the runs. The controller lifts the states with the dictionary,
    condenses the quadratic tracking cost over the horizon with the matrices of the surrogate, solves the bounded least
    squares problem with `scipy.optimize.lsq_linear`, and sets the input to the first value of the solution. The
    plant reads the input through a role, and `CoSimulationSolver` holds it between samples. The plants are those of
    figure 6.6 of the review of Brunton et al.: the forced Duffing oscillator
    `x'' = a x + b x^3 + d x' + f cos(w t) + (3 + cos(x)) u`, the Van der Pol oscillator
    `x'' = -x - mu (x^2 - 1) x' + u`, and the bilinear DC motor, whose second variable is the output that tracks the
    reference. The motor reads four times the input, as Korda and Mezić scale it. One time unit of the system is one
    sample of the controller.

    Args:
        plant:
            `'duffing'`, `'van_der_pol'`, or `'motor'`.

        controlled:
            Whether the controller solves the tracking problem, or holds the input at 0 with no fit.

        dictionary:
            A scikit-learn transformer that computes the observables and keeps the variables among them, or `None` for
            the monomials up to degree 3.

        x:
            The initial value of the first variable.

        y:
            The initial value of the second variable.

        n_runs:
            The number of training runs.

        duration:
            The number of samples of every training run.

        random_state:
            The seed or generator of the training runs.

    Returns:
        system:
            The system, labelled by the plant, with the variables `plant/state.x` and `plant/state.y`, the input
            `controller/input.u`, and the observable `reference`.

    Raises:
        ValidationError:
            The plant is unknown.

    Examples:
        >>> from skcomplex.systems.engineering.control import koopman_model_predictive_control
        >>> system = koopman_model_predictive_control('motor', controlled=False)
        >>> [subsystem.label for subsystem in system.systems]
        ['controller', 'plant']
    """
    if plant not in PLANTS:
        error_msg = f'The plant {plant!r} is unknown. Use one of {", ".join(PLANTS)}.'
        raise ValidationError(error_msg)
    case = PLANTS[plant]
    problem = _fit_problem(case, dictionary, n_runs, duration, random_state) if controlled else None
    controller = MapInteraction(
        'policy',
        target='controller/input',
        state='u',
        sources=['plant/state'],
        function=functools.partial(_hold, problem=problem),
    )
    return System(
        plant,
        systems=[
            System('controller', elements=[ControlInputElement('input')], time=DiscreteTime()),
            System('plant', elements=[PlantStateElement('state', x=x, y=y)], time=ContinuousTime()),
        ],
        interactions=[controller, *_plant_interactions(case, 'plant/state', 'controller/input')],
    )

lqr_oscillator(x=1.0, *, omega=1.0, zeta=0.0, q=1.0, r=0.1)

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.

Parameters:

Name	Type	Description	Default
x	float	The initial displacement.	1.0
omega	float	The natural frequency.	1.0
zeta	float	The damping ratio without control.	0.0
q	float	The weight of the displacement and the velocity in the cost.	1.0
r	float	The weight of the force in the cost.	0.1

Returns:

Name	Type	Description
system	System	The system, with the displacement oscillator.x, the force oscillator.force, and the observable energy.

Examples:

>>> from skcomplex.systems.engineering.control import lqr_oscillator
>>> [round(value, 3) for value in lqr_oscillator().interactions[0].gain]
[2.317, 3.825]

Source code in src/skcomplex/systems/engineering/control/_regulation.py

@register_system(
    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 system, with the displacement `oscillator.x`, the force `oscillator.force`, and the observable
            `energy`.

    Examples:
        >>> from skcomplex.systems.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)
    system = LinearizationAnalysis(['oscillator.force'], state_weights=q * np.eye(2), control_weights=[[r]])
    system.prepare(System(elements=[oscillator])).run()
    gain = system.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 system.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])

minimum_time_control(x=3.0, v=2.0, *, plant='free_particle', interval=0.1)

Return a linear plant steered to the origin of its phase space in minimum time by a sampled bang-bang control.

The plant follows x'' = a x + u with |u| <= 1: a free particle, x'' = u, a harmonic oscillator, x'' + x = u, or an unstable saddle, x'' - x = u. Pontryagin's maximum principle makes the minimum-time control bang-bang, u = ±1, with the switch on the curve of the two trajectories under u = 1 and u = -1 that reach the origin: the control is -1 above the curve and 1 below it. As in Dilão's examples, the controller samples the position and the velocity every interval and holds the control until the next sample, so the plant reaches a neighbourhood of the origin, whose size shrinks with the interval, and stays there. A free particle from (x, v) above the curve reaches the origin at the time v + 2 sqrt(x + v^2 / 2), about 6.47 for the defaults. Only the states between the lines y = -(x + 1) and y = -(x - 1) of the saddle are controllable, and from the others the plant escapes to infinity.

Parameters:

Name	Type	Description	Default
x	float	The initial position.	3.0
v	float	The initial velocity.	2.0
plant	str	The plant, 'free_particle', 'harmonic_oscillator', or 'saddle'.	'free_particle'
interval	float	The time between two samples of the controller.	0.1

Returns:

Name	Type	Description
system	System	The system, with the position plant.x, the velocity plant.y, and the control plant.u, which the samples sample set.

Raises:

Type	Description
ParameterError	The plant is unknown, or the interval is not positive.

Examples:

>>> from skcomplex.systems.engineering.control import minimum_time_control
>>> minimum_time_control(plant='saddle', x=0.5, v=0.0).elements[0].u
-1.0

Source code in src/skcomplex/systems/engineering/control/_minimum_time_control.py

@register_system(
    field='engineering/control',
    tags=('control', 'bang-bang', 'optimal control', 'pontryagin', 'events', 'ode'),
    references=('Pontryagin et al. 1962', 'Dilão 2023'),
    defaults={'until': 10.0, 'dt': 0.01},
)
def minimum_time_control(
    x: float = 3.0,
    v: float = 2.0,
    *,
    plant: str = 'free_particle',
    interval: float = 0.1,
) -> System:
    """Return a linear plant steered to the origin of its phase space in minimum time by a sampled bang-bang control.

    The plant follows `x'' = a x + u` with `|u| <= 1`: a free particle, `x'' = u`, a harmonic oscillator,
    `x'' + x = u`, or an unstable saddle, `x'' - x = u`. Pontryagin's maximum principle makes the minimum-time
    control bang-bang, `u = ±1`, with the switch on the curve of the two trajectories under `u = 1` and `u = -1` that
    reach the origin: the control is -1 above the curve and 1 below it. As in Dilão's examples, the controller samples
    the position and the velocity every `interval` and holds the control until the next sample, so the plant reaches
    a neighbourhood of the origin, whose size shrinks with the interval, and stays there. A free particle from `(x, v)`
    above the curve reaches the origin at the time `v + 2 sqrt(x + v^2 / 2)`, about 6.47 for the defaults. Only the
    states between the lines `y = -(x + 1)` and `y = -(x - 1)` of the saddle are controllable, and from the others the
    plant escapes to infinity.

    Args:
        x:
            The initial position.

        v:
            The initial velocity.

        plant:
            The plant, `'free_particle'`, `'harmonic_oscillator'`, or `'saddle'`.

        interval:
            The time between two samples of the controller.

    Returns:
        system:
            The system, with the position `plant.x`, the velocity `plant.y`, and the control `plant.u`, which the
            samples `sample` set.

    Raises:
        ParameterError:
            The plant is unknown, or the interval is not positive.

    Examples:
        >>> from skcomplex.systems.engineering.control import minimum_time_control
        >>> minimum_time_control(plant='saddle', x=0.5, v=0.0).elements[0].u
        -1.0
    """
    if plant not in PLANTS:
        error_msg = f'Unknown plant {plant!r}. Use one of {", ".join(PLANTS)}.'
        raise ParameterError(error_msg)
    if not interval > 0:
        error_msg = f'The interval of the samples is {interval!r}. Give a positive number.'
        raise ParameterError(error_msg)
    start = float(BangBangPlantElement._control(plant, x, v))
    sample = {'u': lambda s: BangBangPlantElement._control(plant, s.target.x, s.target.y)}
    return System(
        elements=[BangBangPlantElement('plant', x=x, y=v, u=start, a=PLANTS[plant])],
        interactions=[ImpulseInteraction('sample', target='plant', reset=sample, period=interval, start=interval)],
    )

pendulum_swing_up(theta=0.0, omega0=0.5, *, omega=1.5652475842498528, mu=0.1)

Return a pendulum swung up to its inverted position by a bang-bang acceleration of its pivot.

The pivot of an acrobat's pendulum moves horizontally with acceleration f'' = u, so the pendulum follows θ'' + ω^2 sin θ + u cos θ = 0. The control u = -μ sign(θ' cos θ) sign(E_s - E) pumps the energy E = θ'^2 / 2 - ω^2 cos θ towards E_s = ω^2, the energy of the inverted position, from below or from above.

Parameters:

Name	Type	Description	Default
theta	float	The initial angle from the hanging position.	0.0
omega0	float	The initial angular velocity.	0.5
omega	float	The natural frequency sqrt(g / l), by default for g = 9.8 and l = 4.	1.5652475842498528
mu	float	The size of the acceleration of the pivot.	0.1

Returns:

Name	Type	Description
system	System	The system, with the angle pendulum.theta, the position of the pivot pendulum.f, and their velocities.

Examples:

>>> from skcomplex.systems.engineering.control import pendulum_swing_up
>>> pendulum_swing_up(mu=0.2).elements[0].mu
0.2

Source code in src/skcomplex/systems/engineering/control/_control.py

@register_system(
    field='engineering/control',
    tags=('control', 'bang-bang', 'pendulum', 'energy', 'ode'),
    references=('Åström and Furuta 2000', 'Dilão 2023'),
    defaults={'until': 100.0, 'dt': 0.1, 'solver': RK4Solver(dt=0.01)},
)
def pendulum_swing_up(
    theta: float = 0.0,
    omega0: float = 0.5,
    *,
    omega: float = 1.5652475842498528,
    mu: float = 0.1,
) -> System:
    """Return a pendulum swung up to its inverted position by a bang-bang acceleration of its pivot.

    The pivot of an acrobat's pendulum moves horizontally with acceleration `f'' = u`, so the pendulum follows
    `θ'' + ω^2 sin θ + u cos θ = 0`. The control `u = -μ sign(θ' cos θ) sign(E_s - E)` pumps the energy
    `E = θ'^2 / 2 - ω^2 cos θ` towards `E_s = ω^2`, the energy of the inverted position, from below or from above.

    Args:
        theta:
            The initial angle from the hanging position.

        omega0:
            The initial angular velocity.

        omega:
            The natural frequency `sqrt(g / l)`, by default for `g = 9.8` and `l = 4`.

        mu:
            The size of the acceleration of the pivot.

    Returns:
        system:
            The system, with the angle `pendulum.theta`, the position of the pivot `pendulum.f`, and their velocities.

    Examples:
        >>> from skcomplex.systems.engineering.control import pendulum_swing_up
        >>> pendulum_swing_up(mu=0.2).elements[0].mu
        0.2
    """
    pendulum = SwingingPendulumElement('pendulum', theta=theta, theta_d1=omega0, f=0.0, f_d1=0.0, omega=omega, mu=mu)
    return System(elements=[pendulum])

pid_inverted_pendulum(theta=0.1, *, kp=40.0, ki=1.0, kd=8.0, M=1.0, m=0.1, l=0.5)

Return an inverted pendulum on a cart that a PID controller of its angle keeps upright.

The controller sets the force on the cart to kp e + ki ∫ e dt - kd theta' from the error e = -theta of the angle, within the bounds of the force. Near upright the closed loop is stable when kp > (M + m) g and kd > 0, and the angle returns to upright while the cart drifts.

Parameters:

Name	Type	Description	Default
theta	float	The initial angle from upright.	0.1
kp	float	The proportional gain.	40.0
ki	float	The integral gain.	1.0
kd	float	The derivative gain.	8.0
M	float	The mass of the cart.	1.0
m	float	The mass at the end of the rod.	0.1
l	float	The length of the rod.	0.5

Returns:

Name	Type	Description
system	System	The system, with the angle cart.theta, the position cart.x, the force cart.force, and the observable tilt.

Examples:

>>> from skcomplex.systems.engineering.control import pid_inverted_pendulum
>>> pid_inverted_pendulum(kp=50.0).interactions[0].kp
50.0

Source code in src/skcomplex/systems/engineering/control/_regulation.py

@register_system(
    field='engineering/control',
    tags=('control', 'feedback', 'pid', 'unstable equilibrium', 'ode'),
    references=('Åström and Murray 2008',),
    defaults={'until': 5.0, 'dt': 0.05, 'solver': RK4Solver(dt=0.001)},
    observables=lambda system: [ObservableSpec('tilt', _tilt)],
)
def pid_inverted_pendulum(
    theta: float = 0.1,
    *,
    kp: float = 40.0,
    ki: float = 1.0,
    kd: float = 8.0,
    M: float = 1.0,
    m: float = 0.1,
    l: float = 0.5,
) -> System:
    """Return an inverted pendulum on a cart that a PID controller of its angle keeps upright.

    The controller sets the force on the cart to `kp e + ki ∫ e dt - kd theta'` from the error `e = -theta` of the
    angle, within the bounds of the force. Near upright the closed loop is stable when `kp > (M + m) g` and `kd > 0`,
    and the angle returns to upright while the cart drifts.

    Args:
        theta:
            The initial angle from upright.

        kp:
            The proportional gain.

        ki:
            The integral gain.

        kd:
            The derivative gain.

        M:
            The mass of the cart.

        m:
            The mass at the end of the rod.

        l:
            The length of the rod.

    Returns:
        system:
            The system, with the angle `cart.theta`, the position `cart.x`, the force `cart.force`, and the observable
            `tilt`.

    Examples:
        >>> from skcomplex.systems.engineering.control import pid_inverted_pendulum
        >>> pid_inverted_pendulum(kp=50.0).interactions[0].kp
        50.0
    """
    cart = CartPendulumElement('cart', theta=theta, M=M, m=m, l=l)
    controller = PIDInteraction(
        'pid',
        plant='cart',
        integrator='memory',
        control='force',
        measure='theta',
        kp=kp,
        ki=ki,
        kd=kd,
        limits=(-FORCE_LIMIT, FORCE_LIMIT),
    )
    return System(elements=[cart, IntegratorElement('memory')], interactions=[controller])

pid_loop(y=0.0, *, tau=1.0, gain=1.0, Kp=2.0, Ki=1.0, Kd=0.1, setpoint=1.0, change=0.0, disturbance=-0.5)

Return a work process whose output a PID controller keeps at a set point against a disturbance.

The output follows the first-order process τ y' = -y + gain u + d(t), and the controller sets u = Kp e + Ki ∫ e dt - Kd y' from the error e = r(t) - y, with the derivative taken on the output, so (τ + gain Kd) y' = -y + gain (Kp e + Ki ∫ e dt) + d(t). The disturbance d steps in at time 10, and a coordinator moves the set point r by change at time 20. The integral action removes the error that the disturbance leaves under proportional control alone.

Parameters:

Name	Type	Description	Default
y	float	The initial output.	0.0
tau	float	The time constant of the process.	1.0
gain	float	The gain of the actuator.	1.0
Kp	float	The proportional gain.	2.0
Ki	float	The integral gain.	1.0
Kd	float	The derivative gain.	0.1
setpoint	float	The initial set point.	1.0
change	float	The change of the set point at time 20.	0.0
disturbance	float	The size of the disturbance from time 10.	-0.5

Returns:

Name	Type	Description
system	System	The system, with the output process.y and the integral of the error process.integral.

Examples:

>>> from skcomplex.systems.engineering.control import pid_loop
>>> pid_loop(Ki=0.0).elements[0].Ki
0.0

Source code in src/skcomplex/systems/engineering/control/_control.py

@register_system(
    field='engineering/control',
    tags=('control', 'feedback', 'pid', 'governance', 'ode'),
    references=('Mobus 2022', 'Åström and Murray 2008'),
    defaults={'until': 40.0, 'dt': 0.1},
)
def pid_loop(
    y: float = 0.0,
    *,
    tau: float = 1.0,
    gain: float = 1.0,
    Kp: float = 2.0,
    Ki: float = 1.0,
    Kd: float = 0.1,
    setpoint: float = 1.0,
    change: float = 0.0,
    disturbance: float = -0.5,
) -> System:
    """Return a work process whose output a PID controller keeps at a set point against a disturbance.

    The output follows the first-order process `τ y' = -y + gain u + d(t)`, and the controller sets
    `u = Kp e + Ki ∫ e dt - Kd y'` from the error `e = r(t) - y`, with the derivative taken on the output, so
    `(τ + gain Kd) y' = -y + gain (Kp e + Ki ∫ e dt) + d(t)`. The disturbance `d` steps in at time 10, and a
    coordinator moves the set point `r` by `change` at time 20. The integral action removes the error that the
    disturbance leaves under proportional control alone.

    Args:
        y:
            The initial output.

        tau:
            The time constant of the process.

        gain:
            The gain of the actuator.

        Kp:
            The proportional gain.

        Ki:
            The integral gain.

        Kd:
            The derivative gain.

        setpoint:
            The initial set point.

        change:
            The change of the set point at time 20.

        disturbance:
            The size of the disturbance from time 10.

    Returns:
        system:
            The system, with the output `process.y` and the integral of the error `process.integral`.

    Examples:
        >>> from skcomplex.systems.engineering.control import pid_loop
        >>> pid_loop(Ki=0.0).elements[0].Ki
        0.0
    """
    process = ControlledProcessElement(
        'process',
        y=y,
        integral=0.0,
        tau=tau,
        gain=gain,
        Kp=Kp,
        Ki=Ki,
        Kd=Kd,
        setpoint=setpoint,
        change=change,
        disturbance=disturbance,
    )
    return System(elements=[process])

satellite_transfer(s=0.314, velocity=0.0, theta=0.0, *, Lz=0.7, epsilon=0.1, H_final=-0.2)

Return a satellite moved between Kepler orbits of equal angular momentum by radial thrust.

In normalised polar coordinates the satellite follows s'' = L_z^2 / s^3 - 1 / s^2 - ε sigma and θ' = L_z / s^2. Radial thrust keeps the angular momentum and changes the effective energy H = s'^2 / 2 + L_z^2 / (2 s^2) - 1 / s at the rate -ε sigma s', so the controller pushes inwards, sigma = 1, while the satellite falls and outwards, sigma = -1, while it rises, until H reaches the final energy. From the periapsis of the orbit with H = -0.7, the transfer to H = -0.2 takes about 9.81.

Parameters:

Name	Type	Description	Default
s	float	The initial distance to the planet.	0.314
velocity	float	The initial radial velocity.	0.0
theta	float	The initial angle.	0.0
Lz	float	The angular momentum.	0.7
epsilon	float	The thrust.	0.1
H_final	float	The effective energy of the final orbit.	-0.2

Returns:

Name	Type	Description
system	System	The system, with the distance satellite.s, the radial velocity satellite.s_d1, the angle satellite.theta, and the observable energy.

Examples:

>>> from skcomplex.systems.engineering.control import satellite_transfer
>>> satellite_transfer().elements[0].Lz
0.7

Source code in src/skcomplex/systems/engineering/control/_control.py

@register_system(
    field='engineering/control',
    tags=('control', 'orbits', 'celestial mechanics', 'ode'),
    references=('Dilão et al. 2022', 'Dilão 2023'),
    defaults={'until': 20.0, 'dt': 0.05, 'solver': RK4Solver(dt=0.005)},
    observables=lambda system: [ObservableSpec('energy', _energy_of_satellite)],
)
def satellite_transfer(
    s: float = 0.314,
    velocity: float = 0.0,
    theta: float = 0.0,
    *,
    Lz: float = 0.7,
    epsilon: float = 0.1,
    H_final: float = -0.2,
) -> System:
    """Return a satellite moved between Kepler orbits of equal angular momentum by radial thrust.

    In normalised polar coordinates the satellite follows `s'' = L_z^2 / s^3 - 1 / s^2 - ε sigma` and
    `θ' = L_z / s^2`. Radial thrust keeps the angular momentum and changes the effective energy
    `H = s'^2 / 2 + L_z^2 / (2 s^2) - 1 / s` at the rate `-ε sigma s'`, so the controller pushes inwards, `sigma = 1`,
    while the satellite falls and outwards, `sigma = -1`, while it rises, until `H` reaches the final energy. From the
    periapsis of the orbit with `H = -0.7`, the transfer to `H = -0.2` takes about 9.81.

    Args:
        s:
            The initial distance to the planet.

        velocity:
            The initial radial velocity.

        theta:
            The initial angle.

        Lz:
            The angular momentum.

        epsilon:
            The thrust.

        H_final:
            The effective energy of the final orbit.

    Returns:
        system:
            The system, with the distance `satellite.s`, the radial velocity `satellite.s_d1`, the angle
            `satellite.theta`, and the observable `energy`.

    Examples:
        >>> from skcomplex.systems.engineering.control import satellite_transfer
        >>> satellite_transfer().elements[0].Lz
        0.7
    """
    satellite = SatelliteElement(
        'satellite',
        s=s,
        s_d1=velocity,
        theta=theta,
        Lz=Lz,
        epsilon=epsilon,
        H_final=H_final,
    )
    return System(elements=[satellite])

thermostat(T=10.0, *, setting=18.0, furnace=5.0, leak=0.1, outside=10.0, swing=0.0)

Return a room heated by a furnace under a thermostat while heat leaks to the outside.

The furnace adds the gap between the setting and the room temperature, at most furnace degrees per hour, and nothing when the room is warmer than the setting. The room loses leak of its difference to the outside temperature per hour, so it settles below the setting, at (setting + leak * outside) / (1 + leak) for a constant outside. The outside temperature swings by swing around outside over a day, warmest at noon.

Parameters:

Name	Type	Description	Default
T	float	The initial room temperature.	10.0
setting	float	The thermostat setting.	18.0
furnace	float	The largest heating rate of the furnace, in degrees per hour.	5.0
leak	float	The fraction of the difference to the outside lost per hour.	0.1
outside	float	The mean outside temperature.	10.0
swing	float	The amplitude of the daily swing of the outside temperature.	0.0

Returns:

Name	Type	Description
system	System	The system, with the room temperature room.T in hours.

Examples:

>>> from skcomplex.systems.engineering.control import thermostat
>>> thermostat(leak=0.3).elements[0].leak
0.3

Source code in src/skcomplex/systems/engineering/control/_control.py

@register_system(
    field='engineering/control',
    tags=('control', 'feedback', 'stocks and flows', 'ode'),
    references=('Meadows 2008',),
    defaults={'until': 24.0, 'dt': 0.25},
)
def thermostat(
    T: float = 10.0,
    *,
    setting: float = 18.0,
    furnace: float = 5.0,
    leak: float = 0.1,
    outside: float = 10.0,
    swing: float = 0.0,
) -> System:
    """Return a room heated by a furnace under a thermostat while heat leaks to the outside.

    The furnace adds the gap between the setting and the room temperature, at most `furnace` degrees per hour, and
    nothing when the room is warmer than the setting. The room loses `leak` of its difference to the outside
    temperature per hour, so it settles below the setting, at `(setting + leak * outside) / (1 + leak)` for a
    constant outside. The outside temperature swings by `swing` around `outside` over a day, warmest at noon.

    Args:
        T:
            The initial room temperature.

        setting:
            The thermostat setting.

        furnace:
            The largest heating rate of the furnace, in degrees per hour.

        leak:
            The fraction of the difference to the outside lost per hour.

        outside:
            The mean outside temperature.

        swing:
            The amplitude of the daily swing of the outside temperature.

    Returns:
        system:
            The system, with the room temperature `room.T` in hours.

    Examples:
        >>> from skcomplex.systems.engineering.control import thermostat
        >>> thermostat(leak=0.3).elements[0].leak
        0.3
    """
    room = RoomElement('room', T=T, setting=setting, furnace=furnace, leak=leak, outside=outside, swing=swing)
    return System(elements=[room])