Real engineering systems rarely live in a single domain. A servo drive couples electrical storage and dissipation to rotational inertia and elastic transmission. A power system couples electromagnetic energy exchange to mechanical rotor dynamics. A robot manipulator couples motor-side electrical dynamics, rigid-body kinetic energy, elastic joints, and dissipation.

The central idea is to model the system using a stored-energy function \( H(x) \) and a power-conserving interconnection structure. Instead of starting from ad hoc equations, we start from a power balance:

\[ \frac{d}{dt}H(x(t)) = P_{\text{in} }(t) - P_{\text{diss} }(t) \]

This perspective gives physically meaningful states, automatically supports passivity arguments, and makes interconnection of subsystems systematic.

flowchart TD
  A["Identify storage elements in each domain"] --> B["Choose energy states x"]
  B --> C["Write H(x): total stored energy"]
  C --> D["Define effort / flow variables"]
  D --> E["Build power-conserving interconnection"]
  E --> F["Add dissipation and external ports"]
  F --> G["Obtain xdot and output y"]
  G --> H["Check energy balance and passivity"]
  H --> I["Simulate and validate"]
        

Energy-based modeling uses conjugate variables: \( e \) (effort) and \( f \) (flow), whose product is power.

\[ p(t) = e(t)\,f(t) \]

Examples used in this lesson:

• Electrical: effort = voltage \( v \), flow = current \( i \)
• Rotational mechanical: effort = torque \( T \), flow = angular speed \( \omega \)
• Translational mechanical: effort = force \( F \), flow = velocity \( \dot{x} \)
• Hydraulic (review link): effort = pressure, flow = volumetric flow rate
• Thermal (lumped approximation): effort = temperature (or thermodynamic potential), flow = entropy flow / heat-flow proxy

For a vector of domain variables, the instantaneous supplied power is

\[ P = e^\top f \]

If \( x \) collects energy states (fluxes, momenta, generalized coordinates for springs/capacitors, etc.), then the energy gradient \( \nabla H(x) \) generates the conjugate variables. For example:

\[ i = \frac{\partial H}{\partial \phi}, \qquad \omega = \frac{\partial H}{\partial p}, \qquad F_s = \frac{\partial H}{\partial q} \]

where \( \phi \) is magnetic flux linkage, \( p \) is mechanical momentum, and \( q \) is an elastic deformation coordinate.

A standard multi-domain port-Hamiltonian model is

\[ \dot{x} = \big(J(x)-R(x)\big)\nabla H(x) + G(x)u, \qquad y = G(x)^\top \nabla H(x) \]

with structural properties

\[ J(x) = -J(x)^\top, \qquad R(x)=R(x)^\top \succeq 0 \]

Here \( J(x) \) encodes power-conserving interconnection (skew-symmetry ensures zero net internal power), while \( R(x) \) models dissipation.

Proposition (Passivity): The input-output pair \( (u,y) \) is passive with storage function \( H \).

Proof: Differentiate \( H \) along trajectories:

\[ \dot{H} = \nabla H^\top \dot{x} = \nabla H^\top (J-R)\nabla H + \nabla H^\top G u \]

Since \( J=-J^\top \), the quadratic form with \( J \) vanishes:

\[ \nabla H^\top J \nabla H = 0 \]

Therefore

\[ \dot{H} = -\nabla H^\top R \nabla H + \nabla H^\top G u = -\nabla H^\top R \nabla H + y^\top u \le y^\top u \]

because \( R \succeq 0 \). This is exactly the passivity inequality. \(\square\)

Consider a lumped electromechanical actuator with:

• Armature inductance \( L \) and resistance \( R_a \)
• Rotor inertia \( J \) and viscous damping \( b \)
• Torsional shaft stiffness \( k \)
• Electromechanical coupling coefficient \( K_t \)

Choose states \( x = [\phi,\; p,\; q]^\top \), where \( \phi \) is armature flux linkage, \( p \) is angular momentum, and \( q \) is shaft twist.

The stored energy (Hamiltonian) is

\[ H(\phi,p,q)=\frac{\phi^2}{2L}+\frac{p^2}{2J}+\frac{k}{2}q^2 \]

Hence the conjugate variables are

\[ \nabla H= \begin{bmatrix} \phi/L\\ p/J\\ kq \end{bmatrix} = \begin{bmatrix} i\\ \omega\\ T_s \end{bmatrix} \]

A valid port-Hamiltonian realization is

\[ \dot{x} = \left( \underbrace{\begin{bmatrix} 0 & K_t & 0\\ -K_t & 0 & -1\\ 0 & 1 & 0 \end{bmatrix} }_{J_m} - \underbrace{\begin{bmatrix} R_a & 0 & 0\\ 0 & b & 0\\ 0 & 0 & 0 \end{bmatrix} }_{R_m} \right)\nabla H + \underbrace{\begin{bmatrix} 1\\0\\0 \end{bmatrix} }_{G}v \]

where \( v \) is the applied armature voltage. The port output is

\[ y = G^\top \nabla H = i \]

so the input power is \( yv = iv \).

Energy balance: Using the passivity proof above,

\[ \dot{H} = -R_a i^2 - b\omega^2 + iv \]

which states that the rate of stored energy equals electrical input power minus ohmic and viscous losses.

Expanded state equations:

\[ \dot{\phi} = K_t\omega - R_a i + v \]

\[ \dot{p} = -K_t i - kq - b\omega \]

\[ \dot{q} = \omega \]

with \( i=\phi/L \) and \( \omega=p/J \). This model is ideal for simulation and for passivity-based control design.

We now connect the same methodology to a classical power-system dynamic model. For an introductory but rigorous case, consider a single-machine equivalent swing model:

\[ \dot{\delta} = \omega, \qquad M\dot{\omega} = P_m - D\omega - P_{\max}\sin\delta \]

where \( \delta \) is rotor angle deviation, \( \omega \) is speed deviation, \( M \) is inertia, \( D \) is damping, and \( P_m \) is mechanical input power.

Introduce momentum \( p = M\omega \). Define the energy function

\[ H(\delta,p)=\frac{p^2}{2M} + U(\delta), \qquad U(\delta)=P_{\max}(1-\cos\delta)-P_m\delta \]

Then

\[ \frac{\partial H}{\partial p}=\frac{p}{M}=\omega, \qquad \frac{\partial H}{\partial \delta}=P_{\max}\sin\delta - P_m \]

and the dynamics become

\[ \begin{bmatrix} \dot{\delta}\\ \dot{p} \end{bmatrix} = \left( \begin{bmatrix} 0 & 1\\ -1 & 0 \end{bmatrix} - \begin{bmatrix} 0 & 0\\ 0 & D \end{bmatrix} \right) \begin{bmatrix} \partial H/\partial \delta\\ \partial H/\partial p \end{bmatrix} \]

This is a dissipative Hamiltonian form with no external port shown explicitly. If disturbances are modeled as an external power injection \( u \) in the momentum channel, then

\[ \dot{p} = -\frac{\partial H}{\partial \delta} - D\frac{\partial H}{\partial p} + u \]

and the output can be chosen as \( y=\partial H/\partial p=\omega \), so \( yu \) is incremental power.

Local stability insight from energy curvature: At an equilibrium \( (\delta^\star,0) \), we need

\[ P_{\max}\sin\delta^\star = P_m \]

The Hessian of \( H \) is positive definite locally if

\[ \frac{\partial^2 H}{\partial \delta^2}\Big|_{\delta^\star} = P_{\max}\cos\delta^\star > 0 \]

i.e., for stable operating points with \( |\delta^\star| < \pi/2 \). This is an energy-based interpretation of rotor-angle stability near equilibrium.

For an \( n \)-DOF rigid robot manipulator, the Euler-Lagrange equations (with viscous joint dissipation) are

\[ M(q)\ddot{q} + C(q,\dot{q})\dot{q} + D(q)\dot{q} + \nabla V(q) = B(q)u \]

where \( q\in\mathbb{R}^n \) is the joint-coordinate vector, \( M(q) \) is the positive-definite inertia matrix, \( C(q,\dot{q}) \) contains Coriolis/centrifugal effects, \( D(q)\succeq 0 \) is dissipation, and \( V(q) \) is potential energy.

Define momentum coordinates

\[ p = M(q)\dot{q} \]

and Hamiltonian

\[ H(q,p)=\frac{1}{2}p^\top M(q)^{-1}p + V(q) \]

Then the canonical energy relations are

\[ \dot{q} = \frac{\partial H}{\partial p}, \qquad \dot{p} = -\frac{\partial H}{\partial q} - D(q)\frac{\partial H}{\partial p} + B(q)u + \Gamma(q,p) \]

where \( \Gamma(q,p) \) is the skew-power term needed to represent the geometry-dependent Coriolis structure in momentum coordinates. In full port-Hamiltonian form, this is absorbed into a state-dependent skew-symmetric matrix \( J(x) \).

Passivity from actuator torque to joint velocity: If the output is

\[ y = B(q)^\top \frac{\partial H}{\partial p} = B(q)^\top \dot{q} \]

then the robot satisfies

\[ \dot{H} = -\dot{q}^\top D(q)\dot{q} + y^\top u \]

This result is the mathematical backbone behind many passivity-based and energy-shaping controllers used in robotics.

Remark on coordinate consistency: In robotics, using physically meaningful energies is especially important when coupling flexible joints, motor inertias, and links. The multi-domain model often contains motor electrical states, joint elastic states, and rigid-body kinetic energy in one composite Hamiltonian:

\[ H_{\text{total} } = H_{\text{electrical} } + H_{\text{motor} } + H_{\text{joint\;elastic} } + H_{\text{link\;kinetic} } + V_{\text{gravity} } \]

A major advantage of energy-based modeling is modular interconnection. Suppose subsystems \( \Sigma_1, \Sigma_2, \dots, \Sigma_m \) have individual storage functions \( H_i \) and port pairs \( (u_i,y_i) \). If we interconnect internal ports with a skew-symmetric law, the internal power cancels automatically.

Stack all internal ports:

\[ u_{\text{int} } = -K\,y_{\text{int} }, \qquad K=-K^\top \]

Then the internal power is

\[ y_{\text{int} }^\top u_{\text{int} } = -y_{\text{int} }^\top K y_{\text{int} } = 0 \]

so the total storage satisfies

\[ \frac{d}{dt}\sum_{i=1}^{m} H_i \le y_{\text{ext} }^\top u_{\text{ext} } \]

This is the rigorous reason bond-graph/port-Hamiltonian interconnections scale well to complex mechatronic and robotic systems.

flowchart TD
  S1["Subsystem 1 \n(energy storage + \ndissipation)"] --> IC["Power-conserving interconnection"]
  S2["Subsystem 2 \n(energy storage + \ndissipation)"] --> IC
  S3["Subsystem 3 \n(energy storage + \ndissipation)"] --> IC
  IC --> EXT["External ports (actuation / disturbance / measurement)"]
  EXT --> BAL["Total energy balance check"]
  BAL --> SIM["Simulation and controller design"]
        

The following script simulates the mechatronic port-Hamiltonian model using RK4 and numerically verifies the energy balance \( H(t_f)-H(0)\approx \int_0^{t_f}(P_{\text{in} }-P_{\text{diss} })\,dt \).

Chapter18_Lesson5.py

# Chapter18_Lesson5.py
# Port-Hamiltonian mechatronic model (electrical + rotational + elastic)
import numpy as np
import matplotlib.pyplot as plt

L, R_a = 0.35, 1.8
J, b = 0.02, 0.08
k, Kt = 4.0, 0.22

Jm = np.array([[0.0, Kt, 0.0], [-Kt, 0.0, -1.0], [0.0, 1.0, 0.0]])
Rm = np.diag([R_a, b, 0.0])
G = np.array([1.0, 0.0, 0.0])

def gradH(x):
    phi, p, q = x
    return np.array([phi/L, p/J, k*q])

def H(x):
    phi, p, q = x
    return 0.5*phi*phi/L + 0.5*p*p/J + 0.5*k*q*q

def u(t):
    return 8.0 if t < 0.8 else (3.0 if t < 1.6 else 0.0)

def f(t, x):
    return (Jm - Rm) @ gradH(x) + G * u(t)

def rk4(t, x, h):
    k1 = f(t, x)
    k2 = f(t + 0.5*h, x + 0.5*h*k1)
    k3 = f(t + 0.5*h, x + 0.5*h*k2)
    k4 = f(t + h, x + h*k3)
    return x + (h/6.0)*(k1 + 2*k2 + 2*k3 + k4)

t0, tf, h = 0.0, 4.0, 1e-3
t = np.arange(t0, tf + h, h)
x = np.zeros((len(t), 3))
Hhist = np.zeros(len(t)); Pin = np.zeros(len(t)); Pdiss = np.zeros(len(t))
for n in range(len(t)-1):
    g = gradH(x[n]); i, w = g[0], g[1]
    Hhist[n] = H(x[n]); Pin[n] = u(t[n])*i; Pdiss[n] = R_a*i*i + b*w*w
    x[n+1] = rk4(t[n], x[n], h)
g = gradH(x[-1]); Hhist[-1] = H(x[-1]); Pin[-1] = u(t[-1])*g[0]; Pdiss[-1] = R_a*g[0]**2 + b*g[1]**2

lhs = Hhist[-1] - Hhist[0]
rhs = np.trapz(Pin - Pdiss, t)
print("Energy residual =", lhs - rhs)

np.savetxt("Chapter18_Lesson5_python_results.csv",
           np.column_stack([t, x, Hhist, Pin, Pdiss]),
           delimiter=",",
           header="t,phi,p,q,H,Pin,Pdiss", comments="")

plt.figure()
plt.plot(t, x[:,0]/L, label="i(t)")
plt.plot(t, x[:,1]/J, label="omega(t)")
plt.plot(t, x[:,2], label="q(t)")
plt.grid(True); plt.legend(); plt.xlabel("t (s)")
plt.title("Port-Hamiltonian Mechatronic States")
plt.show()
      

Library coverage: numpy (numerics) and matplotlib (plots). This script is suitable for extension to multi-port cases.

The C++ version uses a lightweight RK4 implementation and writes a CSV for later plotting or comparison with Python/MATLAB.

Chapter18_Lesson5.cpp

// Chapter18_Lesson5.cpp
#include <array>
#include <fstream>
#include <iomanip>
#include <iostream>
using Vec3 = std::array<double,3>;

struct P { double L=0.35, Ra=1.8, J=0.02, b=0.08, k=4.0, Kt=0.22; };

double u(double t){ return (t < 0.8) ? 8.0 : ((t < 1.6) ? 3.0 : 0.0); }

Vec3 gradH(const Vec3& x, const P& p){ return {x[0]/p.L, x[1]/p.J, p.k*x[2]}; }

double H(const Vec3& x, const P& p){
    return 0.5*x[0]*x[0]/p.L + 0.5*x[1]*x[1]/p.J + 0.5*p.k*x[2]*x[2];
}

Vec3 f(double t, const Vec3& x, const P& p){
    auto g = gradH(x,p);
    double i=g[0], w=g[1], ks=g[2];
    return {p.Kt*w - p.Ra*i + u(t), -p.Kt*i - ks - p.b*w, w};
}

Vec3 add(const Vec3&a,const Vec3&b,double s){ return {a[0]+s*b[0], a[1]+s*b[1], a[2]+s*b[2]}; }

Vec3 rk4(double t, const Vec3& x, double h, const P& p){
    Vec3 k1=f(t,x,p), k2=f(t+0.5*h, add(x,k1,0.5*h),p), k3=f(t+0.5*h, add(x,k2,0.5*h),p), k4=f(t+h, add(x,k3,h),p);
    return {x[0]+h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6.0,
            x[1]+h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6.0,
            x[2]+h*(k1[2]+2*k2[2]+2*k3[2]+k4[2])/6.0};
}

int main(){
    P p; double h=1e-3, tf=4.0; int N=int(tf/h)+1; Vec3 x{0,0,0};
    std::ofstream out("Chapter18_Lesson5_cpp_results.csv");
    out << "t,phi,p,q,H,Pin,Pdiss\n" << std::setprecision(10);
    double intPow=0.0, prev=0.0; bool first=true; double H0=H(x,p);
    for(int n=0; n < N; ++n){
        double t=n*h; auto g=gradH(x,p); double i=g[0], w=g[1];
        double Pin=u(t)*i, Pdiss=p.Ra*i*i + p.b*w*w, Hx=H(x,p), s=Pin-Pdiss;
        if(!first) intPow += 0.5*h*(prev+s); first=false; prev=s;
        out << t << "," << x[0] << "," << x[1] << "," << x[2] << "," << Hx << "," << Pin << "," << Pdiss << "\n";
        if(n < N-1) x = rk4(t,x,h,p);
    }
    std::cout << "Energy residual = " << (H(x,p)-H0-intPow) << "\n";
    return 0;
}
      

The C++ structure is intentionally minimal and can be embedded in real-time or hardware-in-the-loop workflows.

This Java version follows the same energy-state formulation and produces a CSV file. It is useful for teaching cross-platform simulation and object- oriented modeling of dynamic systems.

Chapter18_Lesson5.java

// Chapter18_Lesson5.java
import java.io.FileWriter;
import java.io.PrintWriter;

public class Chapter18_Lesson5 {
    static class P { double L=0.35, Ra=1.8, J=0.02, b=0.08, k=4.0, Kt=0.22; }

    static double u(double t){ return (t < 0.8) ? 8.0 : ((t < 1.6) ? 3.0 : 0.0); }

    static double[] gradH(double[] x, P p){ return new double[]{x[0]/p.L, x[1]/p.J, p.k*x[2]}; }

    static double H(double[] x, P p){ return 0.5*x[0]*x[0]/p.L + 0.5*x[1]*x[1]/p.J + 0.5*p.k*x[2]*x[2]; }

    static double[] f(double t, double[] x, P p){
        double[] g = gradH(x,p); double i=g[0], w=g[1], ks=g[2];
        return new double[]{p.Kt*w - p.Ra*i + u(t), -p.Kt*i - ks - p.b*w, w};
    }

    static double[] add(double[] a, double[] b, double s){
        return new double[]{a[0]+s*b[0], a[1]+s*b[1], a[2]+s*b[2]};
    }

    static double[] rk4(double t, double[] x, double h, P p){
        double[] k1=f(t,x,p), k2=f(t+0.5*h, add(x,k1,0.5*h),p), k3=f(t+0.5*h, add(x,k2,0.5*h),p), k4=f(t+h, add(x,k3,h),p);
        return new double[]{
            x[0]+h*(k1[0]+2*k2[0]+2*k3[0]+k4[0])/6.0,
            x[1]+h*(k1[1]+2*k2[1]+2*k3[1]+k4[1])/6.0,
            x[2]+h*(k1[2]+2*k2[2]+2*k3[2]+k4[2])/6.0
        };
    }

    public static void main(String[] args) throws Exception {
        P p = new P(); double h=1e-3, tf=4.0; int N=(int)(tf/h)+1; double[] x={0,0,0};
        double intPow=0.0, prev=0.0, H0=H(x,p); boolean first=true;
        try (PrintWriter out = new PrintWriter(new FileWriter("Chapter18_Lesson5_java_results.csv"))) {
            out.println("t,phi,p,q,H,Pin,Pdiss");
            for(int n=0; n < N; n++){
                double t=n*h; double[] g=gradH(x,p); double i=g[0], w=g[1];
                double Pin=u(t)*i, Pdiss=p.Ra*i*i + p.b*w*w, Hx=H(x,p), s=Pin-Pdiss;
                if(!first) intPow += 0.5*h*(prev+s); first=false; prev=s;
                out.printf("%.6f,%.10f,%.10f,%.10f,%.10f,%.10f,%.10f%n", t, x[0], x[1], x[2], Hx, Pin, Pdiss);
                if(n < N-1) x = rk4(t,x,h,p);
            }
        }
        System.out.println("Energy residual = " + (H(x,p)-H0-intPow));
    }
}
      

The Java and C++ implementations are numerically comparable, which helps students validate cross-language consistency of the model.

The MATLAB script uses ode45 and computes the same energy residual. A Simulink realization can be built directly from the equation \( \dot{x}=(J-R)\nabla H + Gu \) using Integrator, Gain, Sum, and Product blocks.

Chapter18_Lesson5.m

% Chapter18_Lesson5.m
% Port-Hamiltonian mechatronic model (electrical + rotational + elastic)
clear; clc;

P.L = 0.35; P.Ra = 1.8; P.J = 0.02; P.b = 0.08; P.k = 4.0; P.Kt = 0.22;
x0 = [0;0;0];
[t,x] = ode45(@(t,x) fph(t,x,P), [0 4], x0);

phi = x(:,1); p = x(:,2); q = x(:,3);
i = phi / P.L; omega = p / P.J;
H = 0.5*phi.^2/P.L + 0.5*p.^2/P.J + 0.5*P.k*q.^2;
u = arrayfun(@uin, t);
Pin = u .* i;
Pdiss = P.Ra*i.^2 + P.b*omega.^2;

lhs = H(end)-H(1);
rhs = trapz(t, Pin-Pdiss);
fprintf('Energy residual = %.3e\n', lhs-rhs);

T = table(t,phi,p,q,H,Pin,Pdiss);
writetable(T,'Chapter18_Lesson5_matlab_results.csv');

plot(t, i, t, omega, t, q); grid on;
legend('i(t)','omega(t)','q(t)'); title('Port-Hamiltonian States');

% Simulink note: implement xdot = (J-R)gradH + G*u using Gain/Sum/Integrator blocks.

function dx = fph(t,x,P)
g = [x(1)/P.L; x(2)/P.J; P.k*x(3)];
dx = [P.Kt*g(2) - P.Ra*g(1) + uin(t);
     -P.Kt*g(1) - g(3) - P.b*g(2);
      g(2)];
end

function v = uin(t)
if t < 0.8
    v = 8.0;
elseif t < 1.6
    v = 3.0;
else
    v = 0.0;
end
end
      

For Simulink, a good exercise is to construct the state derivative with a matrix-based block implementation and then compare against the scripted solution above.

Mathematica is especially useful here because it can both verify the passivity identity symbolically and solve the nonlinear state equations numerically.

Chapter18_Lesson5.nb

Notebook[{
 Cell["Chapter18_Lesson5.nb","Title"],
 Cell["Port-Hamiltonian electromechanical model (electrical + rotational + elastic)","Text"],
 Cell[BoxData@ToBoxes@TraditionalForm[
   H[\[Phi]_,p_,q_] == \[Phi]^2/(2 L) + p^2/(2 J) + (k q^2)/2
 ],"Input"],
 Cell["Symbolic passivity identity","Section"],
 Cell[BoxData@ToBoxes[
   Module[{g,Jm,Rm,G,xdot,Hdot},
    g={phi/L,p/J,k q}; Jm={ {0,Kt,0},{-Kt,0,-1},{0,1,0} }; Rm=DiagonalMatrix[{Ra,b,0}]; G={1,0,0};
    xdot=(Jm-Rm).g + G u;
    Hdot = Expand[g.xdot] // Simplify;
    {"dHdt" -> Hdot, "Target" -> (u (phi/L) - Ra (phi/L)^2 - b (p/J)^2)}
   ]
 ],"Input"],
 Cell["Numerical simulation","Section"],
 Cell[BoxData@ToBoxes[
   Module[{L=0.35,Ra=1.8,J=0.02,b=0.08,k=4.0,Kt=0.22,u,sol,i,w,Hf,pin,pd,lhs,rhs},
    u[t_] := Piecewise[{ {8.0,t<0.8},{3.0,t<1.6} },0.0];
    sol = NDSolveValue[
      {phi'[t]==Kt (p[t]/J)-Ra (phi[t]/L)+u[t],
       p'[t]==-Kt (phi[t]/L)-k q[t]-b (p[t]/J),
       q'[t]==p[t]/J,
       phi[0]==0,p[0]==0,q[0]==0},
      {phi,p,q},{t,0,4}];
    i[t_] := sol[[1]][t]/L; w[t_] := sol[[2]][t]/J;
    Hf[t_] := sol[[1]][t]^2/(2 L) + sol[[2]][t]^2/(2 J) + k sol[[3]][t]^2/2;
    pin[t_] := u[t] i[t]; pd[t_] := Ra i[t]^2 + b w[t]^2;
    lhs = N[Hf[4]-Hf[0]]; rhs = NIntegrate[pin[t]-pd[t],{t,0,4}];
    {"EnergyResidual" -> (lhs-rhs), Plot[{i[t],w[t],sol[[3]][t]},{t,0,4}]}
   ]
 ],"Input"]
}]
      

This notebook format demonstrates symbolic verification and numerical simulation in a single environment.

Suppose we interconnect a motor electrical subsystem, a drivetrain elastic subsystem, and a robot-link mechanical subsystem. A common modeling choice is a composite state \( x = [x_e^\top, x_d^\top, x_r^\top]^\top \) with total storage

\[ H(x)=H_e(x_e)+H_d(x_d)+H_r(x_r) \]

If each subsystem is passive and internal interconnections are power-conserving, the total system remains passive. This gives a direct route to Lyapunov candidates and energy-aware controller design.

In matrix form, a block-structured interconnection is often written as

\[ \dot{x} = \big(J_{\text{blk} }(x)-R_{\text{blk} }(x)\big)\nabla H(x) + G_{\text{blk} }(x)u \]

with \( J_{\text{blk} }=-J_{\text{blk} }^\top \) and \( R_{\text{blk} }\succeq 0 \). The skew-symmetric block entries describe domain couplings (electrical-mechanical, mechanical-elastic, etc.) and guarantee no artificial energy creation.

Problem 1 (Passivity of a port-Hamiltonian model): Let \( \dot{x}=(J-R)\nabla H + Gu \) and \( y=G^\top \nabla H \), with \( J=-J^\top \) and \( R=R^\top \succeq 0 \). Prove that the model is passive from input \( u \) to output \( y \).

Solution: Differentiate the storage:

\[ \dot{H}=\nabla H^\top \dot{x} =\nabla H^\top (J-R)\nabla H + \nabla H^\top G u \]

The skew-symmetric term vanishes because

\[ \nabla H^\top J \nabla H = 0 \]

Hence

\[ \dot{H} = -\nabla H^\top R \nabla H + y^\top u \le y^\top u \]

Therefore the system is passive with storage function \( H \).

Problem 2 (Energy balance for the mechatronic actuator): For the model in Section 4, show explicitly that \( \dot{H} = -R_a i^2 - b\omega^2 + iv \).

Solution: Write

\[ H=\frac{\phi^2}{2L}+\frac{p^2}{2J}+\frac{k}{2}q^2, \qquad \nabla H= \begin{bmatrix} i\\\omega\\kq \end{bmatrix} \]

Using the given pH dynamics,

\[ \dot{H}=\nabla H^\top (J_m-R_m)\nabla H + \nabla H^\top G v \]

Since \( J_m=-J_m^\top \), the interconnection term is zero:

\[ \nabla H^\top J_m \nabla H = 0 \]

Thus

\[ \dot{H}=-\nabla H^\top R_m \nabla H + iv =-(R_a i^2 + b\omega^2) + iv \]

which is the required balance law.

Problem 3 (Swing-equation equilibrium and local minimum of energy): For \( H(\delta,p)=\frac{p^2}{2M}+P_{\max}(1-\cos\delta)-P_m\delta \), determine the equilibrium and the condition for local positive curvature in \( \delta \).

Solution: Equilibrium requires

\[ \dot{\delta}=0 \Rightarrow p^\star = 0 \]

\[ \dot{p}=0 \Rightarrow \frac{\partial H}{\partial \delta}(\delta^\star,p^\star)=0 \Rightarrow P_{\max}\sin\delta^\star = P_m \]

The second derivative is

\[ \frac{\partial^2 H}{\partial \delta^2}=P_{\max}\cos\delta \]

So local convexity (and hence an energy minimum in the angle direction) requires

\[ P_{\max}\cos\delta^\star > 0 \]

For positive \( P_{\max} \), this corresponds to \( |\delta^\star| < \pi/2 \).

Problem 4 (Manipulator passivity with viscous dissipation): For a robot with Hamiltonian \( H(q,p)=\frac{1}{2}p^\top M(q)^{-1}p+V(q) \), assume the damping term is \( D(q)\dot{q} \) with \( D(q)\succeq 0 \). Show that the input torque to output velocity map is passive.

Solution: Let the actuation be \( B(q)u \) and output \( y=B(q)^\top \dot{q} \). In port-Hamiltonian coordinates:

\[ \dot{H} = -\dot{q}^\top D(q)\dot{q} + y^\top u \]

Because \( D(q)\succeq 0 \), the dissipation term is nonnegative:

\[ \dot{q}^\top D(q)\dot{q} \ge 0 \]

Hence

\[ \dot{H} \le y^\top u \]

proving passivity from actuator torque input to joint-velocity output.

Problem 5 (Power-conserving interconnection): Two passive subsystems have internal ports \( (u_1,y_1) \) and \( (u_2,y_2) \). Connect them by

\[ \begin{bmatrix} u_1\\u_2 \end{bmatrix} = \begin{bmatrix} 0 & -1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} y_1\\y_2 \end{bmatrix} \]

Show that the internal interconnection does not create or dissipate power.

Solution: Internal power is

\[ y_1u_1 + y_2u_2 = \begin{bmatrix} y_1 & y_2 \end{bmatrix} \begin{bmatrix} 0 & -1\\ 1 & 0 \end{bmatrix} \begin{bmatrix} y_1\\y_2 \end{bmatrix} \]

The middle matrix is skew-symmetric, so the quadratic form is zero:

\[ y^\top K y = 0 \quad \text{for } K=-K^\top \]

Therefore the interconnection is power-conserving. Any total energy decay is due only to subsystem dissipation, not the coupling itself.

Multi-domain energy-based modeling provides a common mathematical language for mechatronics, power systems, and robotics. The port-Hamiltonian form separates (i) power-conserving interconnection, (ii) dissipation, and (iii) external power exchange, making physical consistency transparent. This structure is not just elegant: it supports stability proofs, modular modeling, and controller design through passivity and energy shaping.