Skip to content

Control

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

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.models.engineering.control import CartPendulumElement
>>> CartPendulumElement('cart', theta=0.1).theta
0.1
Source code in src/skcomplex/models/engineering/control/_cart_pendulum.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
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

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.models.engineering.control import DrivenOscillatorElement
>>> DrivenOscillatorElement('oscillator', x=1.0).omega
1.0
Source code in src/skcomplex/models/engineering/control/_driven_oscillator.py 
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
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

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.models.engineering.control import IntegratorElement
>>> IntegratorElement('memory').value
0.0
Source code in src/skcomplex/models/engineering/control/_integrator.py 
31
32
33
34
35
36
37
38
39
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.models.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/models/engineering/control/_pid.py 
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
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

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.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)
Source code in src/skcomplex/models/engineering/control/_state_feedback.py 
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
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

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 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]
Source code in src/skcomplex/models/engineering/control/_regulation.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
@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])

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 model, with the angle pendulum.theta, the position of the pivot pendulum.f, and their velocities.

Examples:

>>> from skcomplex.models.engineering.control import pendulum_swing_up
>>> pendulum_swing_up(mu=0.2).elements[0].mu
0.2
Source code in src/skcomplex/models/engineering/control/_control.py 
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
@register_model(
    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 model, with the angle `pendulum.theta`, the position of the pivot `pendulum.f`, and their velocities.

    Examples:
        >>> from skcomplex.models.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 model, with the angle cart.theta, the position cart.x, the force cart.force, and the observable tilt.

Examples:

>>> from skcomplex.models.engineering.control import pid_inverted_pendulum
>>> pid_inverted_pendulum(kp=50.0).interactions[0].kp
50.0
Source code in src/skcomplex/models/engineering/control/_regulation.py 
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
@register_model(
    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 model, with the angle `cart.theta`, the position `cart.x`, the force `cart.force`, and the observable
            `tilt`.

    Examples:
        >>> from skcomplex.models.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 model, with the output process.y and the integral of the error process.integral.

Examples:

>>> from skcomplex.models.engineering.control import pid_loop
>>> pid_loop(Ki=0.0).elements[0].Ki
0.0
Source code in src/skcomplex/models/engineering/control/_control.py 
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
@register_model(
    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 model, with the output `process.y` and the integral of the error `process.integral`.

    Examples:
        >>> from skcomplex.models.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 model, with the distance satellite.s, the radial velocity satellite.s_d1, the angle satellite.theta, and the observable energy.

Examples:

>>> from skcomplex.models.engineering.control import satellite_transfer
>>> satellite_transfer().elements[0].Lz
0.7
Source code in src/skcomplex/models/engineering/control/_control.py 
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
@register_model(
    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 model, with the distance `satellite.s`, the radial velocity `satellite.s_d1`, the angle
            `satellite.theta`, and the observable `energy`.

    Examples:
        >>> from skcomplex.models.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 model, with the room temperature room.T in hours.

Examples:

>>> from skcomplex.models.engineering.control import thermostat
>>> thermostat(leak=0.3).elements[0].leak
0.3
Source code in src/skcomplex/models/engineering/control/_control.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
@register_model(
    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 model, with the room temperature `room.T` in hours.

    Examples:
        >>> from skcomplex.models.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])