hydrasim · Electro-Hydraulic Servo Plant Simulation

Electro-Hydraulic Servo Plant Simulation for Control Algorithm Development
Single-joint & Excavator-Arm Models · Validated Physics · Numpy-Only Core

hydrasim is a hydraulic plant simulation package for the HydraServo project, serving as a test bench for control algorithms (ADRC, RLS, self-calibration, global coordination). A decoupled controllers module (PID + LUT velocity feedforward + calibration) is included. Two plant layers:

• Single joint HydraulicPlant — 4-way proportional valve + double-acting asymmetric cylinder, 5-state, validated by 12 physics self-checks.
• Excavator arm ExcavatorArm — 2D rigid-body mechanism (boom/arm/bucket + swing) + 4 valve-controlled actuators + shared variable-displacement pump with LUDV flow distribution, 21-state coupled system, three pluggable integration strategies.

Design principle: the simulation core (numpy-only) is fully decoupled from visualization. The core produces data; viz plugs in via callbacks or recorded data — same pattern as Gym/MuJoCo's step() + optional viewer, so the core can be extracted as a standalone package.

---

1. Installation

Recommended: conda environment.

conda create -n hydrasim python=3.11 -y
conda activate hydrasim
pip install numpy matplotlib dearpygui pytest
pip install -e ".[viz,dev]"

Or use the bundled environment.yml:

conda env create -f environment.yml
conda activate hydrasim
pip install -e .

Verify:

python -m pytest -q          # 90 tests pass
python tools/check_arm_expectations.py   # 9-phase physics expectation check

The simulation core depends only on numpy. matplotlib/dearpygui are only required for the visualization extras.

---

2. Quick Start

Single joint

from hydrasim import HydraulicPlant, HydraulicParams

plant = HydraulicPlant(HydraulicParams())
dt = 1e-3  # 1 kHz control loop
for k in range(1000):
    obs = plant.step(u=0.5, dt=dt)   # u ∈ [-1, 1]
print(obs["x"], obs["P_a"], obs["P_b"])

Excavator arm (default imex stepper)

import numpy as np
from hydrasim import ExcavatorArm, ExcavatorParams

arm = ExcavatorArm(ExcavatorParams(), stepper="imex")
dt = 1e-3
for k in range(1000):
    # u_vec = [boom, arm, bucket, swing], each ∈ [-1, 1]
    obs = arm.step(np.array([0.4, 0.0, 0.0, 0.0]), dt)
print(obs["q"], obs["P_s"], obs["ludv_ratio"])

Core interface (single joint and excavator arm share the same shape):

plant.reset(state=None)              # reset
plant.step(u, dt, F_load=None)        # advance one control step
plant.observe()                       # ground-truth observation
plant.measure(noise=True)             # noisy sensor reading
plant.register_callback(fn)           # broadcast obs after each step

End-effector payload

step(u_vec, dt, tau_load=...) accepts an external joint torque vector, so you can simulate a payload in the bucket:

m_payload = 1000.0   # kg of material in the bucket
r_eff = arm.p.mech.L_boom + arm.p.mech.L_arm * 0.7  # effective lever arm

for k in range(N):
    q = arm.observe()["q"]
    tau_load = np.zeros(4)
    # gravity torque on boom (resists lift): negative sign since tau_ext adds to rhs
    tau_load[0] = -m_payload * 9.81 * r_eff * np.cos(q[0])
    arm.step(u, dt, tau_load=tau_load)

---

3. Gallery

Excavator-arm dig cycle (boom lower → arm/bucket curl → boom raise → swing → dump)

2D linkage animation — full dig cycle, 7.5 s scripted command sequence

Joint angles, dual-chamber pressures, valve commands, and supply pressure over the cycle

Flow saturation (LUDV + P_s droop, C1 coupling)

Four-axis full command + 40% pump displacement → severe flow competition (LUDV ratio ≤ 0.3)

Manifold pressure droops from 35 → 30 MPa; LUDV scales each actuator's supply-side flow proportionally

3D real-time excavator-arm panel (DearPyGui)

Run with python examples/demo_arm_realtime.py — 4-axis sliders + 3D wireframe (orbital camera) + 4 live plots, with single-step debugging

Single-joint open-loop response

Left: step response (u=0.5); right: sine tracking — the basic 5-state HydraulicPlant validates valve dynamics, friction, leakage, etc.

2D cylinder schematic with body / piston / rod / load block

Three steppers are numerically equivalent (TLM vs IMEX)

Single-cylinder gravity-lock test — TLM and IMEX produce identical trajectories, proving the historical "divergence" was bugs, not the integrator

---

4. Three Integration Strategies (steppers)

The numerical challenge of the excavator arm is that a closed hydraulic cylinder is physically a kinematic constraint (joint locked), but is modeled as an extremely stiff spring (β/V). Different steppers handle this differently; switch via the stepper= argument:

stepper	method	locked chamber	speed	use case
"imex"	semi-implicit (explicit q + implicit Pa/Pb 2×2)	stable	fast	default, general use
"tlm"	transmission-line modeling (q/qd and Pa/Pb coupled in a 4×4 linear system)	stable, numerically equivalent to imex	fast	A/B reference
"imex_constrained"	IMEX + kinematic lock on closed-valve joints	stable	fast	backup (not normally needed; for comparison only)

Key finding: TLM ≈ IMEX ≈ Constrained

A/B comparison through the abstraction layer shows the three steppers are numerically equivalent across all tested scenarios — proving that the "locked-chamber drain / divergence" was never an integrator issue, but four independent code bugs:

1. C_leak too large (historical 8e-12 was 5–50× the real value, drained a chamber in ~7 s) → corrected to 1e-13.
2. chamber_volumes(x) semantics mismatch (received cylinder length instead of piston displacement, making Vb negative) → corrected to pass x_cyl - cyl_x_min.
3. Return-side orifice flow set to zero (IMEX implicit pressure update and TLM 4×4 elimination only computed the supply side) → added _orifice_return_linearized.
4. Soft-limit torque sign flipped (v0.3.1): rhs used - tau_soft instead of + tau_soft, turning the restoring force into an accelerating force — a pure positive feedback that caused bucket divergence near the limit. Fixed.

The imex_constrained stepper treats a closed-valve cylinder as a kinematic constraint q_i = const with an adaptive gravity-balance pressure. After the v0.3.1 soft-limit fix, plain IMEX is stable enough; constrained is kept as a more conservative backup (not normally needed).

TLM single-cylinder prototype

The hydrasim.tlm subpackage provides an isolated 5-state single-cylinder prototype (TLMSingleCylinder + IMEXSingleCylinder baseline) for verifying the TLM math in isolation. 6 acceptance tests pass.

---

5. Package Structure

hydraulic_sim/
├── pyproject.toml
├── environment.yml
├── README.md / CHANGELOG.md / LICENSE
├── src/hydrasim/
│   ├── core/                  # numpy-only
│   │   ├── params.py          # HydraulicParams / HydraulicState / FeatureFlags
│   │   ├── valve.py           # deadband, hysteresis, 4-quadrant orifice, LUDV
│   │   ├── friction.py        # Stribeck + viscous friction
│   │   ├── integrators.py     # RK4 / Euler
│   │   └── plant.py           # HydraulicPlant (single joint)
│   ├── mech/                  # 2D planar mechanism
│   │   ├── linkage.py         # M(q)/C(q,q̇)/G(q)/FK + soft limits + soil force
│   │   ├── joint_map.py       # cylinder↔joint IK / Jacobian (triangle linkage)
│   │   └── motor.py           # valve-controlled hydraulic motor
│   ├── hydraulics/            # pump and supply manifold
│   │   ├── pump.py            # variable-displacement pump (constant power)
│   │   └── manifold.py        # P_s capacitive state + LUDV distribute()
│   ├── steppers/              # pluggable integration strategies
│   │   ├── base.py            # ArmStepper ABC
│   │   ├── _helpers.py        # shared orifice linearization / 2×2 implicit update
│   │   ├── imex.py            # IMEXStepper
│   │   ├── tlm.py             # TLMStepper
│   │   └── constrained.py     # ConstrainedIMEXStepper
│   ├── tlm/                   # TLM single-cylinder prototype
│   │   └── single_cylinder.py
│   ├── params/excavator.py    # MechParams / PumpParams / MotorParams / ExcavatorParams
│   ├── controllers/           # decoupled controller layer (step(obs,dt)->u)
│   │   ├── base.py            # Controller ABC
│   │   ├── pid.py             # PIDController (anti-windup, D-filter)
│   │   ├── feedforward.py     # LUT / analytic / directional velocity FF
│   │   ├── velocity.py        # CylinderVelocityController (FF+PID)
│   │   ├── calibration.py     # velocity LUT calibration
│   │   └── utils.py           # RateLimiter, FirstOrderLowPass
│   ├── arm.py                 # ExcavatorArm top-level coupled object
│   ├── data/recorder.py       # ring buffer + npz save/load
│   └── viz/                   # optional visualization (matplotlib + DearPyGui)
│       ├── replay.py / animator.py / animator_arm.py / realtime.py
│       ├── cylinder_view.py   # 2D single-cylinder schematic
│       └── viewer3d.py        # 3D wireframe excavator view
├── docs/                      # user docs (physics baseline/constants, controllers, benchmarks) + README assets
├── examples/                  # 5 demos
├── tests/                     # 14 test files, 90 tests
└── tools/                     # check_arm_expectations + plot/tune scripts + diag/

---

6. Physical Model

Single joint state (5-D)

[x, v, P_a, P_b, x_v]   # displacement, velocity, rodless-chamber pressure,
                        # rod-side-chamber pressure, valve spool position

Mechanics: m·v̇ = P_a·A_a − P_b·A_b − F_friction(v) − F_load

Chamber continuity: Ṗ_a = β/V_a·(Q_a − Q_leak − A_a·v), Ṗ_b = β/V_b·(Q_b + Q_leak + A_b·v)

4-quadrant orifice: Q = K_v·|x_v|·sign(ΔP)·√(2/ρ·|ΔP|) (x_v≥0 extend: P_s→A, B→P_t)

Hydraulic natural frequency: ω_h = √(β_e/m·(A_a²/V_a + A_b²/V_b)), ~26 Hz at mid-stroke by default.

Excavator arm state (21-D)

[q(4), qd(4), boom(Pa,Pb,xv), arm(...), bucket(...), motor(...), P_s]

Cylinder/motor displacement x and velocity v are derived from IK(q) and J·q̇ (not integrated independently) to avoid algebraic loops.

Two core couplings:

• Mechanism coupling: the planar serial chain M(q) has off-diagonal terms; moving one joint changes the gravity torque on the others.
• Flow coupling (C1): four actuators share one pump; under compound motion they compete for flow, and P_s (supply manifold capacitive state) droops dynamically. Under saturation, LUDV scales each actuator's supply-side flow by r = Q_pump/ΣQ_demand.

Default parameter magnitudes (20-ton excavator, calibrated to PC200-8)

parameter	value	note
bore (boom/arm/bucket)	120/135/115 mm	Arm > Boom > Bucket (PC200-8)
P_s	25 MPa	supply pressure
Q_nom	1.5e-3 m³/s	rated flow (~90 L/min)
β_e	1e8 Pa	engineering value with entrained air (pure oil 1.4e9)
C_leak	1e-13 m³/s/Pa	cross-chamber leakage (realistic magnitude)
tau_v	12 ms	valve spool first-order time constant
link lengths	5.7/2.9/1.5 m	boom/arm/bucket (PC200-8)
I_swing	1.5e5 kg·m²	swing inertia

---

7. Toggleable Nonlinearities

p = HydraulicParams().with_flags(deadband=False, hysteresis=False, friction=True)

flag	meaning
deadband	valve deadband
hysteresis	backlash hysteresis
valve_dynamics	valve spool first-order dynamics
friction	Stribeck + viscous friction
leakage	cross-chamber internal leakage
position_volume	stroke-dependent chamber volume
sensor_noise	sensor noise / quantization

---

8. Data Recording & Visualization

from hydrasim.data import DataRecorder

rec = DataRecorder(capacity=200_000, fields=("t", "q", "P_a", "P_b", "P_s"))
arm.register_callback(rec)
# ... run simulation ...
rec.save("run.npz")

Offline plotting/animation:

from hydrasim.viz import plot_timeseries, plot_arm_timeseries, animate_cylinder, animate_arm
plot_arm_timeseries(data, save="arm.png", show=False)
animate_arm(data, arm.p, save="arm.gif", show=False)

Real-time panels (DearPyGui):

python examples/demo_realtime.py            # single joint, 2D cylinder schematic
python examples/demo_arm_realtime.py        # excavator arm, 3D wireframe + curves

---

9. Demos

# single joint
python examples/demo_open_loop.py step --save        # step response
python examples/demo_open_loop.py sine --save        # sine tracking
python examples/demo_realtime.py                     # DearPyGui real-time panel

# excavator arm
python examples/demo_excavator_arm.py --save          # open-loop dig cycle
python examples/demo_excavator_arm.py --save --saturate  # flow saturation (LUDV + P_s droop)
python examples/demo_arm_realtime.py                 # 3D real-time panel

# TLM vs IMEX comparison
python examples/demo_tlm_vs_imex.py --save           # single-cylinder gravity-lock

Outputs go to examples/_out/*.{png,gif,npz}. tools/check_arm_expectations.py runs the 9-phase physics expectation check (dig 5 + sat 4, all pass).

---

10. Validation

90 passed

test file	#	coverage
test_plant.py	12	single-joint physics self-checks
test_controllers_pid.py	13	PID controller (anti-windup, D-filter)
test_velocity_controller.py	8	velocity inner-loop basics
test_velocity_tracking_performance.py	8	P1 sine / load / reversal / calibration
test_arm.py	7	excavator arm IMEX regression
test_arm_constrained.py	6	constrained stepper
test_scaling_laws.py	6	scaling laws
test_tlm_single.py	6	TLM single-cylinder prototype
test_analytical_benchmarks.py	5	analytical benchmarks
test_conservation_laws.py	5	conservation laws
test_control_axis.py	5	control-axis response
test_swing_brake.py	4	swing brake
test_counterbalance_valve.py	3	counterbalance valve
test_geom_sensitivity.py	2	geometry sensitivity
check_arm_expectations.py	9 phases	dig 5/5 + sat 4/4 physics expectations

---

11. Known Limitations & Future Work

This is an MVP plant for control-algorithm development, not a high-fidelity AMESim-grade hydraulic simulator.

Numerical-method limitations:

1. Coriolis via numerical differentiation: coriolis_gravity calls mass_matrix 6× per step; ~30–60 s wall-clock per 1 s simulated. Increasing substep_dt (2e-4 still stable) helps; analytical ∂M/∂q is a future optimization.
2. Constrained stepper uses quasi-static lock: closed-valve chamber pressure tracks the current pose's gravity balance; it does not model the small elastic leakage of a real hydraulic lock. Sufficient for a control test bench; high-fidelity use needs a DAE constraint solver.
3. 2D planar: swing is decoupled from planar geometry; no 3D Coriolis coupling.

Model limitations: 4. Friction is tanh-smoothed Stribeck, not true stick-slip; can upgrade to Karnopp/LuGre. 5. Hysteresis is a single backlash operator; can upgrade to Prandtl-Ishlinskii for self-calibration. 6. β_e is constant by default; beta_e_fn(state) is reserved for online oil-temperature/air-fraction estimation. 7. Links are lumped-inertia, rigid, no flex/backlash; cylinder↔joint is a 1-DOF triangle mapping, no parallel closed chain. 8. Soil force is a vertical spring/damper at the bucket tip, not real digging resistance. 9. Supply manifold is a single lumped volume; LUDV is ideal proportional distribution, no pressure-compensator dynamics. 10. Bucket cylinder is a direct triangle drive; the real bell-crank + short link amplifier is not yet modeled (quasi-static equivalent is a planned follow-up). 11. Payload is applied as an external torque only; it does not yet augment I_bucket/I_boom (rotational inertia of a loaded bucket).

Upgrade paths:

• For higher-fidelity stiff handling: scipy Radau + analytical Jacobian, or Assimulo/Sundials IDA (DAE constraints).
• For a mature hydraulic component library: OpenHydraulics + OpenModelica (FMI → Python), or Project Chrono (pychrono).

---

License

MIT, see LICENSE.

Changelog

See CHANGELOG.md.