DC motor with an inertia
A DC motor whose armature circuit of components drives an inertia on a Lagrangian coordinate.
- Name:
dc_motor - Field:
engineering/mechatronics - Formalism: ordinary differential equations
- References: Fritzson 2015
Select this system
from skcomplex.systems import create_simulation, create_system, describe
system = create_system('dc_motor')
entry = describe('dc_motor')
trajectory = create_simulation('dc_motor').prepare().run().trajectory_
Factory
@register_system(
field='engineering/mechatronics',
tags=('mechatronics', 'circuit', 'acausal components', 'euler-lagrange', 'ode'),
references=('Fritzson 2015',),
defaults={'until': 10.0, 'dt': 0.01},
observables=lambda system: [
ObservableSpec('speed', _speed_of_rotor),
ObservableSpec('current', _current_of_winding),
],
)
def dc_motor(
V: float = 12.0,
*,
R: float = 1.0,
L: float = 0.5,
k: float = 0.5,
J: float = 0.1,
b: float = 0.2,
) -> System:
"""Return a DC motor whose armature circuit of components drives an inertia on a Lagrangian coordinate.
A voltage source drives a resistor, an inductor, and the armature in series. The armature has the voltage `k phi'`,
and its current exerts the torque `k i` on the angle of an inertia with viscous friction, so the angular velocity
settles at `k V / (R b + k^2)`.
Args:
V:
The voltage of the source.
R:
The resistance of the armature circuit.
L:
The inductance of the armature circuit.
k:
The motor constant.
J:
The moment of inertia of the rotor.
b:
The coefficient of viscous friction of the rotor.
Returns:
system:
The system, with the angle `rotor.phi`, the angular velocity `rotor.phi_d1`, and the current `winding.p_i`.
Examples:
>>> from skcomplex.systems.engineering.mechatronics import dc_motor
>>> dc_motor(k=0.3).interactions[-1].k
0.3
"""
return System(
elements=[
VoltageSourceElement('source', V=V),
ResistorElement('resistor', R=R),
InductorElement('coil', L=L),
ArmatureElement('winding'),
GroundElement('ground'),
InertiaElement('rotor', J=J, b=b),
],
interactions=[
ConnectionInteraction('supply', members=['source', 'resistor'], ports=['p', 'p']),
ConnectionInteraction('series', members=['resistor', 'coil'], ports=['n', 'p']),
ConnectionInteraction('terminal', members=['coil', 'winding'], ports=['n', 'p']),
ConnectionInteraction('neutral', members=['winding', 'source', 'ground'], ports=['n', 'n', 'p']),
EulerLagrangeInteraction(
'rotation',
target='rotor',
states='phi',
lagrangian=_rotor_lagrangian,
dissipation=_rotor_dissipation,
),
MotorCouplingInteraction('motor', armature='winding', inertia='rotor', k=k),
],
)
Components
Elements:
VoltageSourceElement, fromElementResistorElement, fromElementInductorElement, fromElementArmatureElement, fromElementGroundElement, fromElementInertiaElement, fromElement
Interactions:
ConnectionInteractionEulerLagrangeInteractionMotorCouplingInteraction, fromInteraction
class VoltageSourceElement(TwoPinElement):
"""Ideal source of the constant voltage `V = p_v - n_v`.
Args:
label:
The label of the element, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
p_v:
The potential of the positive pin.
p_i:
The current into the positive pin.
n_v:
The potential of the negative pin.
n_i:
The current into the negative pin.
V:
The voltage.
Examples:
>>> from skcomplex.systems.engineering.common.electrical import VoltageSourceElement
>>> VoltageSourceElement('battery', V=9.0).V
9.0
"""
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,
V: Annotated[object, ParameterSpec()] = 1.0,
) -> None:
super().__init__(label, visual, p_v=p_v, p_i=p_i, n_v=n_v, n_i=n_i)
self.V = V
def _get_branch_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the voltage of the source.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The equation law of the branch.
"""
return [EquationLaw(s.p_v - s.n_v, s.V)]
class ResistorElement(TwoPinElement):
"""Resistor that obeys Ohm's law, `p_v - n_v = R p_i`.
Args:
label:
The label of the element, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
p_v:
The potential of the positive pin.
p_i:
The current into the positive pin.
n_v:
The potential of the negative pin.
n_i:
The current into the negative pin.
R:
The resistance.
Examples:
>>> from skcomplex.systems.engineering.common.electrical import ResistorElement
>>> ResistorElement('load', R=2.0).R
2.0
"""
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,
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.R = R
def _get_branch_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return Ohm's law.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The equation law of the branch.
"""
return [EquationLaw(s.p_v - s.n_v, s.R * s.p_i)]
class InductorElement(TwoPinElement):
"""Inductor of current `i`, with `L i' = p_v - n_v` and `i = p_i`.
Args:
label:
The label of the element, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
p_v:
The potential of the positive pin.
p_i:
The current into the positive pin.
n_v:
The potential of the negative pin.
n_i:
The current into the negative pin.
i:
The initial current.
L:
The inductance.
Examples:
>>> from skcomplex.systems.engineering.common.electrical import InductorElement
>>> InductorElement('coil', L=0.5).L
0.5
"""
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,
i: Annotated[object, StateSpec()] = 0.0,
L: 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.i = i
self.L = L
def _get_branch_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the voltage and the current of the inductor.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The equation laws of the branch.
"""
return [EquationLaw(s.L * derivative(s.i, 1), s.p_v - s.n_v), EquationLaw(s.i, s.p_i)]
class ArmatureElement(TwoPinElement):
"""Winding of a motor, a two-pin element whose voltage the coupling with the inertia of the motor gives.
Args:
label:
The label of the element, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
p_v:
The potential of the positive pin.
p_i:
The current into the positive pin.
n_v:
The potential of the negative pin.
n_i:
The current into the negative pin.
Examples:
>>> from skcomplex.systems.engineering.mechatronics import ArmatureElement
>>> sorted(ArmatureElement.ports)
['n', 'p']
"""
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,
) -> None:
super().__init__(label, visual, p_v=p_v, p_i=p_i, n_v=n_v, n_i=n_i)
class GroundElement(Element):
"""Reference of the potentials of a circuit, whose one pin has zero potential.
Args:
label:
The label of the element, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
p_v:
The potential of the pin.
p_i:
The current into the pin.
Examples:
>>> from skcomplex.systems.engineering.common.electrical import GroundElement
>>> sorted(GroundElement.ports)
['p']
"""
ports: ClassVar[Mapping[str, ConnectorSpec]] = {'p': ELECTRICAL}
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
p_v: Annotated[object, StateSpec(algebraic=True)] = 0.0,
p_i: Annotated[object, StateSpec(algebraic=True)] = 0.0,
) -> None:
super().__init__(label, visual=visual)
self.p_v = p_v
self.p_i = p_i
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the zero potential of the pin."""
return [EquationLaw(s.p_v, 0)]
class InertiaElement(Element):
"""Rotating inertia of the angle `phi`, a coordinate of a Lagrangian.
Args:
label:
The label of the element, unique among the parts of its system.
visual:
How the element is drawn, or `None` for `default_visual_`.
phi:
The initial angle.
phi_d1:
The initial angular velocity.
J:
The moment of inertia.
b:
The coefficient of viscous friction.
Examples:
>>> from skcomplex.systems.engineering.mechatronics import InertiaElement
>>> InertiaElement('rotor', J=0.5).J
0.5
"""
def __init__(
self,
label: str,
visual: VisualSpec | None = None,
*,
phi: Annotated[object, StateSpec(order=2)] = 0.0,
phi_d1: object = 0.0,
J: Annotated[object, ParameterSpec(domain=POSITIVE)] = 1.0,
b: Annotated[object, ParameterSpec(domain=POSITIVE)] = 1.0,
) -> None:
super().__init__(label, visual=visual)
self.phi = phi
self.phi_d1 = phi_d1
self.J = J
self.b = b
class ConnectionInteraction(Interaction):
"""Node of connected ports, whose potentials are equal and whose flows sum to zero.
The interaction binds the elements of one node through a variadic role and names the port of each. A node with many
ports is one connection, and a port belongs to one connection. A port that no connection names has a zero flow,
which the compiler adds.
Args:
label:
The label of the interaction, unique among the parts of its system.
members:
The paths of the elements of the node.
ports:
The name of the port of every element, in the same order.
Raises:
ValidationError:
The members and the ports are empty or differ in length.
Examples:
>>> from typing import Annotated, ClassVar
>>> from skcomplex.core.components.elements import Element
>>> from skcomplex.core.components.interactions import ConnectionInteraction
>>> from skcomplex.core.components.laws import EquationLaw
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.specs import ConnectorSpec, ParameterSpec, StateSpec
>>> PIN = ConnectorSpec('pin', potential='v', flow='i')
>>> class ConductanceElement(Element):
... ports: ClassVar = {'a': PIN}
...
... def __init__(
... self,
... label,
... a_v: Annotated[float, StateSpec(algebraic=True)] = 0.0,
... a_i: Annotated[float, StateSpec(algebraic=True)] = 0.0,
... g: Annotated[float, ParameterSpec()] = 1.0,
... e: Annotated[float, ParameterSpec()] = 0.0,
... ):
... super().__init__(label)
... self.a_v = a_v
... self.a_i = a_i
... self.g = g
... self.e = e
...
... def _get_laws(self, s):
... return [EquationLaw(s.a_i, s.g * (s.a_v - s.e))]
>>> system = System(
... elements=[ConductanceElement('first', e=1.0), ConductanceElement('second', e=3.0)],
... interactions=[ConnectionInteraction('node', members=['first', 'second'], ports=['a', 'a'])],
... )
>>> trajectory = Simulation(system, dt=1.0).prepare().run(until=1.0).trajectory_
>>> float(trajectory.values('first.a_v')[-1])
2.0
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (EquationLaw,)
def __init__(
self,
label: str,
*,
members: Annotated[Sequence[str], RoleSpec(variadic=True, elements=(ElementKind.ANY,))],
ports: Sequence[str],
) -> None:
super().__init__(label)
self.members = members
self.ports = ports
def _check(self) -> None:
"""Raise an error for members and ports that are empty or differ in length.
Raises:
ValidationError:
The type declares states, or the members and the ports are empty or differ in length.
"""
super()._check()
if not self.members or len(self.members) != len(self.ports):
error_msg = (
f'{self.label!r} connects {len(self.members)} members through {len(self.ports)} ports. Name one port '
'per member.'
)
raise ValidationError(error_msg)
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the equality of the potentials of the ports and the balance of their flows.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the members.
Returns:
laws:
The equation laws of the node.
"""
potentials, flows = [], []
for member, port in zip(namespace.members, self.ports, strict=True):
connector = member.ports[port]
potentials.append(getattr(member, f'{port}_{connector.potential}'))
flows.append(getattr(member, f'{port}_{connector.flow}'))
equalities: list[BasePart] = [EquationLaw(potentials[0], other) for other in potentials[1:]]
return [*equalities, EquationLaw(sum(flows, SR(0)), 0)]
class EulerLagrangeInteraction(Interaction):
"""Equations of motion `d/dt dL/dq' - dL/dq = Q - dR/dq'` of the coordinates `q` of a Lagrangian `L(ns)`.
The coordinates are states of order 2, and the Lagrangian, the Rayleigh dissipation function `R`, and the
generalized forcing `Q` read them and their first derivatives `derivative(q, 1)`. The interaction returns one
equation law per coordinate, which the compiler solves jointly for the second derivatives, so a Lagrangian with a
mass matrix, such as a double pendulum's, needs no inversion by hand. The equations of one coordinate from several
interactions add, so their Lagrangians, dissipations, and forcings sum: a coupling, such as a spring between two
pendulums, is an interaction of its own whose Lagrangian holds only the coupling energy. The coordinates are the
states of every role that binds one element or group, `target` by default, and every component of a vector state
is a coordinate.
A Lagrangian, dissipation, or forcing function takes the namespace of the interaction, with the namespaces of the
roles, the time `t`, and the symbols of the parameters a subclass declares. The forcing returns one generalized
force per coordinate, in the order of the roles and of the states, or one expression for a single coordinate.
Args:
label:
The label of the interaction, unique among the parts of its system.
target:
The path of the element or group whose states are the coordinates.
states:
The name of the state of order 2 that is a coordinate, or the names of several states.
lagrangian:
A function of the namespace that returns the Lagrangian.
dissipation:
A function of the namespace that returns the Rayleigh dissipation function, or `None` for none.
forcing:
A function of the namespace that returns the generalized forces, or `None` for none.
Raises:
ValidationError:
A state is not a name, or a function is not callable.
Examples:
>>> from typing import Annotated
>>> from sage.all__sagemath_symbolics import cos
>>> from skcomplex.core.components.elements import Element
>>> from skcomplex.core.components.interactions import EulerLagrangeInteraction
>>> from skcomplex.core.components.laws.utils import derivative
>>> from skcomplex.core.components.system import System
>>> from skcomplex.core.engine.runtime import Simulation
>>> from skcomplex.core.specs import StateSpec
>>> class AngleElement(Element):
... def __init__(self, label, q: Annotated[float, StateSpec(order=2)] = 0.1, q_d1: float = 0.0):
... super().__init__(label)
... self.q = q
... self.q_d1 = q_d1
>>> swing = EulerLagrangeInteraction(
... 'swing',
... target='bob',
... states='q',
... lagrangian=lambda s: derivative(s.target.q, 1) ** 2 / 2 + cos(s.target.q),
... )
>>> system = System(elements=[AngleElement('bob')], interactions=[swing])
>>> trajectory = Simulation(system, dt=0.01).prepare().run(until=3.2).trajectory_
>>> round(float(trajectory.values('bob.q')[-1]), 2)
-0.1
"""
law_classes: ClassVar[tuple[type[BaseLaw], ...]] = (EquationLaw,)
def __init__(
self,
label: str,
*,
target: Annotated[str, RoleSpec(elements=(ElementKind.ANY,), states=('states',))],
states: str | Sequence[str] = 'x',
lagrangian: Callable[[SimpleNamespace], object] | None = None,
dissipation: Callable[[SimpleNamespace], object] | None = None,
forcing: Callable[[SimpleNamespace], object] | None = None,
) -> None:
super().__init__(label)
self.target = target
self.states = states
self.lagrangian = lagrangian
self.dissipation = dissipation
self.forcing = forcing
def _get_names(self) -> tuple[str, ...]:
"""Return the names of the states that are coordinates."""
return (self.states,) if isinstance(self.states, str) else tuple(self.states)
def _check(self) -> None:
"""Raise an error for a state that is not a name, or a function that is not callable.
Raises:
ValidationError:
The type declares states, a state is not a name, or a function is not callable.
"""
super()._check()
names = self._get_names()
if not names or not all(isinstance(name, str) and name.isidentifier() for name in names):
error_msg = f'{type(self).__name__} {self.label!r} reads the states {self.states!r}. Give names of states.'
raise ValidationError(error_msg)
for name in ('lagrangian', 'dissipation', 'forcing'):
value = getattr(self, name)
if value is not None and not callable(value):
error_msg = f'The {name} of {self.label!r} is {value!r}. Give a function of the namespace.'
raise ValidationError(error_msg)
def _get_lagrangian(self, namespace: SimpleNamespace) -> object:
"""Return the Lagrangian, the hook a subclass overrides to derive it from its parameters."""
if self.lagrangian is None:
error_msg = f'{type(self).__name__} {self.label!r} has no Lagrangian. Give a function of the namespace.'
raise ValidationError(error_msg)
return self.lagrangian(namespace)
def _get_coordinates(self, namespace: SimpleNamespace) -> list[Expression]:
"""Return the coordinates, every component of the states of every role that binds one path."""
coordinates = []
for name, spec in self._get_specs().items():
if not isinstance(spec, RoleSpec) or spec.variadic or spec.noise:
continue
for state in self._get_names():
value = getattr(getattr(namespace, name), state)
coordinates += list(value) if isinstance(value, FreeModuleElement) else [value]
return coordinates
def _get_forces(self, namespace: SimpleNamespace, count: int) -> list[Expression]:
"""Return the generalized force of every coordinate, zero without a forcing."""
if self.forcing is None:
return [SR(0)] * count
forces = self.forcing(namespace)
forces = list(forces) if isinstance(forces, FreeModuleElement | list | tuple) else [forces]
if len(forces) != count:
error_msg = f'The forcing of {self.label!r} returns {len(forces)} forces for {count} coordinates.'
raise ValidationError(error_msg)
return [SR(force) for force in forces]
def _get_laws(self, namespace: SimpleNamespace) -> list[BasePart]:
"""Return the Euler-Lagrange equation of every coordinate.
Args:
namespace:
The namespace of the symbols the laws can use, with the namespaces of the roles.
Returns:
laws:
The equation laws of the coordinates.
"""
coordinates = self._get_coordinates(namespace)
velocities = {derivative(q, 1): SR.var(f'_velocity_{index}') for index, q in enumerate(coordinates)}
back = {symbol: item for item, symbol in velocities.items()}
slots = list(velocities.values())
lagrangian = SR(self._get_lagrangian(namespace)).subs(velocities)
dissipation = SR(0 if self.dissipation is None else self.dissipation(namespace)).subs(velocities)
forces = self._get_forces(namespace, len(coordinates))
time = SR(namespace.t)
laws: list[BasePart] = []
for q, slot, force in zip(coordinates, slots, forces, strict=True):
momentum = lagrangian.diff(slot)
change = momentum.diff(time) + sum(
momentum.diff(other) * other_slot + momentum.diff(other_slot) * derivative(other, 2)
for other, other_slot in zip(coordinates, slots, strict=True)
)
lhs = (change - lagrangian.diff(q)).subs(back)
rhs = (force - dissipation.diff(slot)).subs(back)
laws.append(EquationLaw(lhs == rhs, target=q))
return laws
class MotorCouplingInteraction(Interaction):
"""Coupling of an armature and an inertia by a motor constant `k`.
The voltage of the armature is the back electromotive force `k phi'`, and its current exerts the torque `k p_i` on
the angle, which adds to the Euler-Lagrange equation of the angle.
Args:
label:
The label of the interaction, unique among the parts of its system.
armature:
The path of the armature.
inertia:
The path of the inertia.
k:
The motor constant.
Examples:
>>> from skcomplex.systems.engineering.mechatronics import MotorCouplingInteraction
>>> MotorCouplingInteraction('motor', armature='winding', inertia='rotor', k=0.5).k
0.5
"""
def __init__(
self,
label: str,
*,
armature: Annotated[str, RoleSpec()],
inertia: Annotated[str, RoleSpec()],
k: Annotated[object, ParameterSpec(domain=POSITIVE)] = 1.0,
) -> None:
super().__init__(label)
self.armature = armature
self.inertia = inertia
self.k = k
def _get_laws(self, s: SimpleNamespace) -> list[BaseLaw]:
"""Return the back electromotive force and the torque on the angle.
Args:
s:
The namespace of the symbols the laws can use.
Returns:
laws:
The equation laws of the coupling.
"""
return [
EquationLaw(s.armature.p_v - s.armature.n_v, s.k * derivative(s.inertia.phi, 1)),
EquationLaw(-s.k * s.armature.p_i, 0, target=s.inertia.phi),
]
Equations
source.p_v: diff(p_v_source(t), t) == 0
source.p_i: diff(p_i_source(t), t) == 0
source.n_v: diff(n_v_source(t), t) == 0
source.n_i: diff(n_i_source(t), t) == 0
resistor.p_v: diff(p_v_resistor(t), t) == 0
resistor.p_i: diff(p_i_resistor(t), t) == 0
resistor.n_v: diff(n_v_resistor(t), t) == 0
resistor.n_i: diff(n_i_resistor(t), t) == 0
coil.p_v: diff(p_v_coil(t), t) == 0
coil.p_i: diff(p_i_coil(t), t) == 0
coil.n_v: diff(n_v_coil(t), t) == 0
coil.n_i: diff(n_i_coil(t), t) == 0
coil.i: diff(i_coil(t), t) == -(R_resistor*i_coil(t) + k_motor*diff(phi_rotor(t), t) - V_source)/L_coil
winding.p_v: diff(p_v_winding(t), t) == 0
winding.p_i: diff(p_i_winding(t), t) == 0
winding.n_v: diff(n_v_winding(t), t) == 0
winding.n_i: diff(n_i_winding(t), t) == 0
ground.p_v: diff(p_v_ground(t), t) == 0
ground.p_i: diff(p_i_ground(t), t) == 0
rotor.phi: diff(phi_rotor(t), t, t) == (k_motor*i_coil(t) - b_rotor*diff(phi_rotor(t), t))/J_rotor