A digital twin is a simulation model that remains synchronized with a physical robot by continuously ingesting real sensor data. Unlike a one-off simulator, a twin is state-aligned and parameter-aligned with the real system. It supports prediction, diagnosis, and rapid “what-if” experimentation.

In this introductory course, we build twins of subsystems (e.g., a single actuator or a mobile base) rather than full kinematic/dynamic models of complex arms.

flowchart TD
  A["Physical robot subsystem"] --> B["Sensors stream y(t)"]
  B --> C["State estimator"]
  C --> D["Digital twin state x_hat(t)"]
  D --> E["Physics simulator"]
  E --> F["Predicted outputs y_hat(t)"]
  F --> G["Residual r(t)=y-y_hat"]
  G --> H["Parameter update"]
  H --> E
  D --> I["Planner / monitor / what-if runs"]
        

The workflow above captures the key loop: measure → estimate state → simulate → compare → update parameters → resimulate.

Let the real robot subsystem be described by a continuous-time model with state \( x(t) \in \mathbb{R}^n \), input \( u(t) \in \mathbb{R}^m \), and output \( y(t) \in \mathbb{R}^p \):

\[ \dot{x}(t) = f(x(t), u(t), \theta), \qquad y(t) = h(x(t), \theta) \]

Here \( \theta \in \mathbb{R}^q \) collects physical parameters (mass, friction, motor constants, etc.). The digital twin maintains an internal estimate \( \hat{x}(t) \) and \( \hat{\theta}(t) \):

\[ \dot{\hat{x}}(t) = f(\hat{x}(t), u(t), \hat{\theta}(t)) \]

Define the state error \( e_x(t)=x(t)-\hat{x}(t) \) and output residual \( r(t)=y(t)-\hat{y}(t) \), where \( \hat{y}(t)=h(\hat{x}(t), \hat{\theta}(t)) \). A fundamental objective of the twin is:

\[ \lim_{t\to\infty} \|e_x(t)\| = 0, \qquad \lim_{t\to\infty} \|r(t)\| = 0 \]

To keep the lesson within linear-control background, we consider a local linearization near an operating point, yielding:

\[ \dot{x}(t) = A x(t) + B u(t), \qquad y(t)=C x(t) \]

The twin runs the same model with estimated state:

\[ \dot{\hat{x}}(t) = A \hat{x}(t) + B u(t) \]

If we inject measurement feedback as a Luenberger-corrector (standard in linear control), the twin becomes:

\[ \dot{\hat{x}}(t) = A \hat{x}(t) + B u(t) + L\big(y(t)-C\hat{x}(t)\big) \]

The error dynamics are then:

\[ \dot{e}_x(t) = (A-LC)e_x(t) \]

Proposition (Twin convergence): If the pair \( (A,C) \) is observable, then there exists a gain \( L \) such that \( A-LC \) is Hurwitz. Hence \( e_x(t)\to 0 \).

Proof sketch: Observability implies controllability of \( (A^\top, C^\top) \). By pole placement on the dual system, choose \( L \) so the eigenvalues of \( A-LC \) lie in the open left-half plane. Then the linear system \( \dot{e}_x=(A-LC)e_x \) is exponentially stable, giving the stated limit. □

Parameter synchronization is performed by minimizing residual energy over a time window \( [t_0, t_0+T] \):

\[ \hat{\theta} = \arg\min_{\theta} \int_{t_0}^{t_0+T} \|y(t)-\hat{y}(t;\theta)\|^2 \, dt \]

In discrete time (sampling period \( \Delta t \)), define stacked residuals \( r_k = y_k - \hat{y}_k \). A basic least-squares update uses:

\[ \hat{\theta} = \arg\min_{\theta}\sum_{k=0}^{N-1}\|r_k(\theta)\|^2 \]

When the output is linear in parameters, i.e., \( \hat{y}_k = \Phi_k \theta \), closed-form LS gives:

\[ \hat{\theta} = \left(\sum_{k}\Phi_k^\top \Phi_k\right)^{-1} \left(\sum_{k}\Phi_k^\top y_k\right) \]

We now describe a small, implementable workflow for a robot subsystem.

1. Choose subsystem and signals. Pick something with accessible sensors, e.g., a single wheel drive or a 1-DOF joint. Identify inputs \( u \) and outputs \( y \).
2. Build physics model. Write \( f(\cdot) \) and \( h(\cdot) \). Use parameters \( \theta \).
3. Implement simulator. Numerically integrate: \( \hat{x}_{k+1}=\hat{x}_k+\Delta t\, f(\hat{x}_k,u_k,\hat{\theta}) \).
4. Synchronize state. Update \( \hat{x} \) using measurement correction (Luenberger/Kalman-style) so residuals shrink.
5. Identify parameters. Periodically solve LS or gradient descent for \( \hat{\theta} \).
6. Validate twin. Compare predicted vs real trajectories; accept if errors are within tolerance.
7. Use twin. Run predictions, fault detection, or sim-to-real tuning.

flowchart TD
  S0["1. Select subsystem"] --> S1["2. Derive model f,h"]
  S1 --> S2["3. Implement simulator"]
  S2 --> S3["4. State sync (observer)"]
  S3 --> S4["5. Parameter sync (LS)"]
  S4 --> S5["6. Validate vs data"]
  S5 --> S6["7. Deploy twin for what-if / monitoring"]
  S6 --> S3
        

Consider a DC motor driving a wheel. The standard linear motor model is:

\[ J\dot{\omega} + b\,\omega = K_t i - \tau_L, \qquad L\dot{i} + R i = V - K_e \omega \]

where \( \omega \) is wheel angular speed, \( i \) armature current, and parameters \( \theta = (J,b,K_t,K_e,R,L) \). The load torque \( \tau_L \) is treated as a disturbance.

Choose state \( x=[\omega\;\; i]^\top \), input \( u=V \), output \( y=\omega \). Then:

\[ \dot{x} = \begin{bmatrix} -\frac{b}{J} & \frac{K_t}{J} \\ -\frac{K_e}{L} & -\frac{R}{L} \end{bmatrix} x + \begin{bmatrix} 0 \\ \frac{1}{L} \end{bmatrix} u + \begin{bmatrix} -\frac{1}{J} \\ 0 \end{bmatrix}\tau_L \]

The simplest twin assumes \( \tau_L=0 \), and relies on state correction from measurements to absorb mismatch.

Suppose we only measure \( \omega \). Then \( C=[1\;\;0] \) and we can design a gain \( L=[\ell_1\;\;\ell_2]^\top \) so that \( A-LC \) is stable (from Section 2).

For parameter identification, note that the first equation can be rearranged:

\[ \dot{\omega} = -\frac{b}{J}\omega + \frac{K_t}{J} i \]

Let \( \theta_1=b/J \) and \( \theta_2=K_t/J \). With numerical derivative \( \dot{\omega}_k \approx (\omega_{k+1}-\omega_k)/\Delta t \), we obtain a linear regression:

\[ \dot{\omega}_k = -\theta_1 \omega_k + \theta_2 i_k \quad\Rightarrow\quad \dot{\omega}_k = \Phi_k \theta,\;\; \Phi_k = [-\omega_k\;\; i_k] \]

Thus the LS formula in Section 2 directly updates \( (\theta_1,\theta_2) \).

5.1 Python (NumPy + Simple Euler Simulator)

import numpy as np

# ---------- True (physical) parameters ----------
J_true, b_true, Kt_true = 0.02, 0.1, 0.05
Ke_true, R_true, L_true = 0.05, 2.0, 0.5

# ---------- Twin parameters (initial guesses) ----------
J_hat, b_hat, Kt_hat = 0.03, 0.2, 0.04
Ke_hat, R_hat, L_hat = Ke_true, R_true, L_true  # assume known for simplicity

dt = 0.01
T = 5.0
N = int(T/dt)

# Input voltage
t = np.arange(N)*dt
V = 6.0 * (t > 0.5)  # step after 0.5s

# Storage
x_true = np.zeros((N,2))  # [omega, i]
x_hat  = np.zeros((N,2))
y_meas = np.zeros(N)

# Observer gain (chosen stable)
L_obs = np.array([30.0, 5.0])  # [l1, l2]^T

def f_motor(x, V, J, b, Kt, Ke, R, L):
    omega, i = x
    domega = (-b/J)*omega + (Kt/J)*i
    di     = (-Ke/L)*omega - (R/L)*i + (1.0/L)*V
    return np.array([domega, di])

# ---------- Simulate physical system and twin ----------
for k in range(N-1):
    # physical model
    x_true[k+1] = x_true[k] + dt*f_motor(x_true[k], V[k],
                                        J_true, b_true, Kt_true, Ke_true, R_true, L_true)
    y_meas[k] = x_true[k,0] + np.random.normal(0, 0.02)  # measure omega

    # twin predictor
    x_hat[k+1] = x_hat[k] + dt*f_motor(x_hat[k], V[k],
                                      J_hat, b_hat, Kt_hat, Ke_hat, R_hat, L_hat)

    # correction (Luenberger)
    residual = y_meas[k] - x_hat[k,0]
    x_hat[k+1] += dt*L_obs*residual

# ---------- Least-squares update for theta1=b/J and theta2=Kt/J ----------
omega = y_meas
i_hat = x_hat[:,1]
domega = np.diff(omega)/dt
Phi = np.column_stack((-omega[:-1], i_hat[:-1]))
theta_ls = np.linalg.inv(Phi.T@Phi) @ (Phi.T@domega)

theta1_hat, theta2_hat = theta_ls
print("Estimated theta1=b/J:", theta1_hat)
print("Estimated theta2=Kt/J:", theta2_hat)

# Update twin parameters using J_hat (keeping J_hat fixed here)
b_hat  = theta1_hat * J_hat
Kt_hat = theta2_hat * J_hat
print("Updated b_hat, Kt_hat:", b_hat, Kt_hat)
      

This script creates a toy digital twin for a motorized wheel. The observer stabilizes state mismatch, and LS reduces parameter mismatch.

5.2 C++ (Eigen + Euler Integration)

#include <iostream>
#include <Eigen/Dense>
#include <random>

using Eigen::Vector2d;
using Eigen::MatrixXd;

Vector2d f_motor(const Vector2d& x, double V,
                 double J, double b, double Kt,
                 double Ke, double R, double L) {
    double omega = x(0), i = x(1);
    double domega = (-b/J)*omega + (Kt/J)*i;
    double di     = (-Ke/L)*omega - (R/L)*i + (1.0/L)*V;
    return Vector2d(domega, di);
}

int main() {
    // True parameters
    double Jt=0.02, bt=0.1, Ktt=0.05, Ket=0.05, Rt=2.0, Lt=0.5;
    // Twin parameters
    double Jh=0.03, bh=0.2, Kth=0.04, Keh=Ket, Rh=Rt, Lh=Lt;

    double dt=0.01, T=5.0;
    int N = (int)(T/dt);

    Vector2d x_true(0,0), x_hat(0,0);
    Vector2d Lobs(30.0, 5.0);

    std::default_random_engine rng(0);
    std::normal_distribution<double> noise(0.0, 0.02);

    MatrixXd Phi(N-1, 2);
    Eigen::VectorXd domega(N-1);

    double omega_prev = 0.0;

    for(int k=0;k<N-1;k++){
        double t = k*dt;
        double V = (t>0.5)?6.0:0.0;

        // Physical step
        x_true += dt * f_motor(x_true, V, Jt, bt, Ktt, Ket, Rt, Lt);
        double y_meas = x_true(0) + noise(rng);

        // Twin predict
        x_hat += dt * f_motor(x_hat, V, Jh, bh, Kth, Keh, Rh, Lh);

        // Correct
        double r = y_meas - x_hat(0);
        x_hat += dt * Lobs * r;

        // Build LS regression for domega = -theta1*omega + theta2*i
        if(k>0){
            domega(k-1) = (y_meas - omega_prev)/dt;
            Phi(k-1,0)  = -omega_prev;
            Phi(k-1,1)  = x_hat(1);
        }
        omega_prev = y_meas;
    }

    // Least squares theta = (Phi^T Phi)^-1 Phi^T domega
    MatrixXd A = Phi.transpose()*Phi;
    Eigen::VectorXd bvec = Phi.transpose()*domega;
    Eigen::VectorXd theta = A.ldlt().solve(bvec);

    std::cout << "theta1=b/J estimate: " << theta(0) << std::endl;
    std::cout << "theta2=Kt/J estimate: " << theta(1) << std::endl;
    return 0;
}
      

5.3 Java (Minimal Matrix Ops)

import java.util.Random;

public class MotorTwin {
    static double[] fMotor(double[] x, double V,
                           double J, double b, double Kt,
                           double Ke, double R, double L){
        double omega = x[0], i = x[1];
        double domega = (-b/J)*omega + (Kt/J)*i;
        double di     = (-Ke/L)*omega - (R/L)*i + (1.0/L)*V;
        return new double[]{domega, di};
    }

    public static void main(String[] args){
        double Jt=0.02, bt=0.1, Ktt=0.05, Ket=0.05, Rt=2.0, Lt=0.5;
        double Jh=0.03, bh=0.2, Kth=0.04, Keh=Ket, Rh=Rt, Lh=Lt;

        double dt=0.01, T=5.0;
        int N=(int)(T/dt);

        double[] xTrue={0,0}, xHat={0,0};
        double[] Lobs={30.0, 5.0};

        Random rng=new Random(0);

        double[][] Phi=new double[N-1][2];
        double[] domega=new double[N-1];

        double omegaPrev=0.0;

        for(int k=0;k<N-1;k++){
            double t=k*dt;
            double V=(t>0.5)?6.0:0.0;

            double[] dxTrue=fMotor(xTrue,V,Jt,bt,Ktt,Ket,Rt,Lt);
            xTrue[0]+=dt*dxTrue[0]; xTrue[1]+=dt*dxTrue[1];

            double yMeas=xTrue[0]+0.02*rng.nextGaussian();

            double[] dxHat=fMotor(xHat,V,Jh,bh,Kth,Keh,Rh,Lh);
            xHat[0]+=dt*dxHat[0]; xHat[1]+=dt*dxHat[1];

            double r=yMeas-xHat[0];
            xHat[0]+=dt*Lobs[0]*r; xHat[1]+=dt*Lobs[1]*r;

            if(k>0){
                domega[k-1]=(yMeas-omegaPrev)/dt;
                Phi[k-1][0]=-omegaPrev;
                Phi[k-1][1]=xHat[1];
            }
            omegaPrev=yMeas;
        }

        // Solve normal equations for 2x2 case
        double a11=0,a12=0,a22=0,b1=0,b2=0;
        for(int k=0;k<N-2;k++){
            double p1=Phi[k][0], p2=Phi[k][1];
            a11+=p1*p1; a12+=p1*p2; a22+=p2*p2;
            b1+=p1*domega[k]; b2+=p2*domega[k];
        }
        double det=a11*a22-a12*a12;
        double theta1=( a22*b1-a12*b2)/det;
        double theta2=(-a12*b1+a11*b2)/det;

        System.out.println("theta1=b/J estimate: "+theta1);
        System.out.println("theta2=Kt/J estimate: "+theta2);
    }
}
      

5.4 MATLAB / Simulink

In MATLAB, we can implement the twin via state-space simulation plus a correction term. The LS update reuses matrix formulas directly.

% True parameters
Jt=0.02; bt=0.1; Ktt=0.05; Ket=0.05; Rt=2.0; Lt=0.5;
% Twin guesses
Jh=0.03; bh=0.2; Kth=0.04; Keh=Ket; Rh=Rt; Lh=Lt;

dt=0.01; T=5; N=T/dt;
t=(0:N-1)*dt;
V=6*(t>0.5);

xTrue=zeros(2,N);
xHat=zeros(2,N);
yMeas=zeros(1,N);
Lobs=[30;5];

% Dynamics function
fMotor=@(x,V,J,b,Kt,Ke,R,L)[(-b/J)*x(1)+(Kt/J)*x(2);
                           (-Ke/L)*x(1)-(R/L)*x(2)+(1/L)*V];

for k=1:N-1
    xTrue(:,k+1)=xTrue(:,k)+dt*fMotor(xTrue(:,k),V(k),Jt,bt,Ktt,Ket,Rt,Lt);
    yMeas(k)=xTrue(1,k)+0.02*randn;

    xHat(:,k+1)=xHat(:,k)+dt*fMotor(xHat(:,k),V(k),Jh,bh,Kth,Keh,Rh,Lh);
    r=yMeas(k)-xHat(1,k);
    xHat(:,k+1)=xHat(:,k+1)+dt*Lobs*r;
end

omega=yMeas;
iHat=xHat(2,:);
domega=diff(omega)/dt;
Phi=[-omega(1:end-1)'  iHat(1:end-1)'];
thetaLS=(Phi'*Phi)\(Phi'*domega');

theta1=thetaLS(1); theta2=thetaLS(2);
disp(['theta1=b/J estimate: ', num2str(theta1)])
disp(['theta2=Kt/J estimate: ', num2str(theta2)])

bh=theta1*Jh; Kth=theta2*Jh; %#ok<NASGU>
      

Simulink sketch: Build two parallel blocks: (1) “Plant” using true parameters; (2) “Twin” using guessed parameters. Feed measured \( \omega \) into a summing junction to form \( r=y-\hat{y} \), then multiply by gains \( \ell_1,\ell_2 \) and add into the twin’s integrators.

Problem 1 (Observer Stability): Consider the linear digital twin \( \dot{\hat{x}} = A\hat{x}+Bu+L(y-C\hat{x}) \). Show that if \( A-LC \) is Hurwitz, then \( \|x(t)-\hat{x}(t)\| \to 0 \).

Solution: Define \( e_x=x-\hat{x} \). Subtracting twin from plant gives \( \dot{e}_x=(A-LC)e_x \). Since \( A-LC \) is Hurwitz, the origin of this LTI system is exponentially stable, hence \( e_x(t)=e^{(A-LC)t}e_x(0)\to 0 \).

Problem 2 (Least Squares Identification): Let measurements satisfy \( \dot{\omega}_k = -\theta_1 \omega_k + \theta_2 i_k \). Derive the normal equations for \( \hat{\theta} \).

Solution: Stack \( N \) samples: \( \mathbf{d} = \Phi \theta \) where \( \mathbf{d}=[\dot{\omega}_0,\dots,\dot{\omega}_{N-1}]^\top \) and \( \Phi=[-\omega_k\;\; i_k]_{k=0}^{N-1} \). The LS cost is \( J(\theta)=\|\mathbf{d}-\Phi\theta\|^2 \). Setting gradient to zero:

\[ \frac{\partial J}{\partial \theta} = -2\Phi^\top(\mathbf{d}-\Phi\theta)=0 \quad\Rightarrow\quad \Phi^\top \Phi \hat{\theta} = \Phi^\top \mathbf{d} \]

giving \( \hat{\theta}=(\Phi^\top\Phi)^{-1}\Phi^\top\mathbf{d} \) whenever \( \Phi^\top\Phi \) is invertible.

Problem 3 (Residual-Based Validation): Assume residuals are evaluated as \( r_k = y_k-\hat{y}_k \). Propose a quantitative acceptance test for the twin.

Solution: A simple test uses normalized RMS error:

\[ \text{NRMSE} = \frac{\sqrt{\frac{1}{N}\sum_{k=0}^{N-1} r_k^2}} {\max_k |y_k| - \min_k |y_k|} \]

Accept the twin if \( \text{NRMSE} < \epsilon \) for a domain-specific tolerance \( \epsilon \) (e.g., \( 0.05 \) for 5% relative error).

Problem 4 (Discrete-Time Twin): For the subsystem \( \dot{x}=Ax+Bu \), derive the forward-Euler twin update.

Solution: Forward Euler approximates \( \dot{x}(t_k)\approx (x_{k+1}-x_k)/\Delta t \), giving:

\[ \hat{x}_{k+1} = \hat{x}_k + \Delta t (A\hat{x}_k + Bu_k) = (I+\Delta t\,A)\hat{x}_k + \Delta t\,B u_k \]

Problem 5 (Twin Sensitivity): Suppose the twin uses incorrect friction \( b \). Show that the predicted steady-state speed under constant voltage differs proportionally to \( 1/b \).

Solution: At steady state, \( \dot{\omega}=0 \) and \( \dot{i}=0 \). From the motor model:

\[ 0 = -b\omega + K_t i \quad\Rightarrow\quad i = \frac{b}{K_t}\omega \]

Substitute into the electrical steady-state equation \( 0 = -K_e \omega - R i + V \):

\[ 0 = -K_e\omega - R\frac{b}{K_t}\omega + V \quad\Rightarrow\quad \omega_\infty = \frac{V}{K_e + (Rb/K_t)} \]

For fixed \( V,K_e,R,K_t \), the term \( Rb/K_t \) scales linearly in \( b \), so \( \omega_\infty \propto 1/b \) for moderate friction-dominant regimes.

We built a simple but complete digital-twin workflow for a robot subsystem. The twin is formalized as a simulator synchronized by (i) state correction via an observer and (ii) parameter correction via least squares. This closes the sim-to-real loop introduced earlier and provides a reusable template for more complex twins in later courses.