Minimum-time control of a free particle, a harmonic oscillator, and a saddle
A linear plant steered to the origin of its phase space in minimum time by a sampled bang-bang control.
- Name:
minimum_time_control - Field:
engineering/control - Formalism: ordinary differential equations
- References: Pontryagin et al. 1962; Dilão 2023
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('minimum_time_control')
entry = describe('minimum_time_control')
trajectory = create_simulation('minimum_time_control').prepare().run().trajectory_
Factory
@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)],
)
Components
Elements:
BangBangPlantElement, fromElement
Interactions:
ImpulseInteraction
class BangBangPlantElement(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.
Args:
label:
The label of the plant.
visual:
How the plant is drawn, or `None` for `default_visual_`.
x:
The initial position.
y:
The initial velocity.
u:
The initial control, -1 or 1.
a:
The coefficient of the position in the acceleration, 0, -1, or 1.
Examples:
>>> from skcomplex.systems.engineering.control import BangBangPlantElement
>>> BangBangPlantElement('plant', a=-1.0).a
-1.0
"""
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
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the motion of the plant under the held control.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The laws of the component.
"""
return [RateLaw(s.x, s.y), RateLaw(s.y, s.a * s.x + s.u), RateLaw(s.u, 0)]
@staticmethod
def _switching_curve(plant: str, x: Expression) -> Expression:
"""Return the velocity `g(x)` of the switching curve of a plant at a position.
Args:
plant:
The plant, `'free_particle'`, `'harmonic_oscillator'`, or `'saddle'`.
x:
The position, a number or a symbol.
Returns:
velocity:
The velocity of the switching curve.
Raises:
ParameterError:
The plant is unknown.
"""
if plant == 'free_particle':
curve: Expression = -sign(x) * sqrt(2 * abs_symbolic(x))
elif plant == 'harmonic_oscillator':
left = sqrt(abs_symbolic(1 - (x + 1) ** 2))
right = -sqrt(abs_symbolic(1 - (x - 1) ** 2))
curve = where(abs_symbolic(x) > CIRCLES_END, 0, where(x < 0, left, right))
elif plant == 'saddle':
curve = where(x <= 0, sqrt(abs_symbolic(x**2 - 2 * x)), -sqrt(abs_symbolic(x**2 + 2 * x)))
else:
error_msg = f'Unknown plant {plant!r}. Use one of {", ".join(PLANTS)}.'
raise ParameterError(error_msg)
return curve
@staticmethod
def _control(plant: str, x: Expression, y: Expression) -> Expression:
"""Return the bang-bang control of a plant, -1 above its switching curve and 1 below it.
Args:
plant:
The plant, `'free_particle'`, `'harmonic_oscillator'`, or `'saddle'`.
x:
The position.
y:
The velocity.
Returns:
control:
The control.
"""
control: Expression = where(y > BangBangPlantElement._switching_curve(plant, x), -1, 1)
return control
class ImpulseInteraction(Interaction):
"""Impulses that set states of the target to new values at given times, or at periodic times.
The impulses fall at every time of `times` and, with a `period`, at `start + j * period` for `j` from 0, `count`
times or without end. An impulse at the initial time of a run does not fire. The new values are functions of the
namespace, which holds `namespace.target`, the states and parameters of the target, and the time `namespace.t`. An
impulse sets the states of every member of a group at once, and the trajectory records every impulse among its
events, under the label of the interaction.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states the impulses set.
reset:
A mapping from the name of a state of the target to the function `f(namespace)` of its new value.
times:
The times of the impulses.
period:
The period of the periodic impulses, or `None` for none.
start:
The time of the first periodic impulse.
count:
The number of periodic impulses, or `None` for impulses without end.
Raises:
ValidationError:
The impulses set no state, a new value is not a function, the impulses have no time, the period is not
positive, or the count is not a positive integer.
Examples:
>>> from skcomplex.core.components.interactions import ImpulseInteraction
>>> reset = {'c': lambda s: s.target.c + 1.0}
>>> dose = ImpulseInteraction('dose', target='body', reset=reset, period=8.0, count=3)
>>> dose.start
0.0
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (CrossingLaw,)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('reset',))],
reset: Mapping[str, Callable[[SimpleNamespace], object]],
times: Sequence[float] = (),
period: float | None = None,
start: float = 0.0,
count: int | None = None,
) -> None:
super().__init__(label)
self.target = target
self.reset = reset
self.times = times
self.period = period
self.start = start
self.count = count
def _check(self) -> None:
"""Raise an error for impulses without states or times, or with a period or a count they cannot take.
Raises:
ValidationError:
The impulses set no state, a new value is not a function, the impulses have no time, the period is not
positive, or the count is not a positive integer.
"""
super()._check()
if not self.reset or not all(callable(function) for function in self.reset.values()):
error_msg = (
f'The impulses {self.label!r} set {self.reset!r}. Map state names to functions of the namespace.'
)
raise ValidationError(error_msg)
if not self.times and self.period is None:
error_msg = f'The impulses {self.label!r} have no time. Give `times`, or a `period`.'
raise ValidationError(error_msg)
if self.period is not None and not self.period > 0:
error_msg = f'The period of the impulses {self.label!r} is {self.period!r}. Give a positive number.'
raise ValidationError(error_msg)
if self.count is not None and (
isinstance(self.count, bool) or not isinstance(self.count, int) or self.count < 1
):
error_msg = f'The count of the impulses {self.label!r} is {self.count!r}. Give a positive integer, or None.'
raise ValidationError(error_msg)
def _get_periodic_guard(self, time: Expression, period: float) -> Expression:
"""Return a guard that crosses zero at every periodic impulse and stays away from zero outside them."""
guard = where(time < self.start - period / 2, -1, sin(pi * (time - self.start) / period))
if self.count is None:
return guard
last = self.count - SR(1) / 2
return where(time > self.start + last * period, sin(pi * last), guard)
def _get_laws(self, namespace: SimpleNamespace) -> list[BaseLaw]:
"""Return the crossing laws of time at the impulses.
Args:
namespace:
The namespace of the symbols the laws can use, with `namespace.target` and the time `namespace.t`.
Returns:
laws:
One crossing law per given time and one for the periodic impulses.
"""
reset = {
getattr(self._get_member(namespace, 'target'), name): SR(f(namespace)) for name, f in self.reset.items()
}
laws: list[BaseLaw] = [CrossingLaw(namespace.t - time, reset, 1, label=self.label) for time in self.times]
if self.period is not None:
laws.append(CrossingLaw(self._get_periodic_guard(namespace.t, self.period), reset, 0, label=self.label))
return laws
Equations
plant.x: diff(x_plant(t), t) == y_plant(t)
plant.y: diff(y_plant(t), t) == a_plant*x_plant(t) + u_plant(t)
plant.u: diff(u_plant(t), t) == 0