Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

4.5 Dynamic Feasibility Constraints

A trajectory must obey the model

Dynamic feasibility requires the proposed state and control histories to satisfy the governing equations and initial conditions:

x˙(t)f(t,x,z,u,p,c,a,d)=0,x(t0)=x0.\dot x(t)-f(t,x,z,u,p,c,a,d)=0,\qquad x(t_0)=x_0.

In simulation-based nested optimization, an integrator enforces this relation implicitly for each trial design. In simultaneous optimization, state samples are decision variables and numerical defect equations enforce the dynamics explicitly.

For forward Euler on a mesh tkt_k,

ζk=xk+1xkΔtkf(tk,xk,uk,p)=0.\zeta_k=x_{k+1}-x_k-\Delta t_k f(t_k,x_k,u_k,p)=0.

Higher-order transcription schemes use more accurate quadrature and interpolation, but the idea is the same: a small objective value is meaningless if the trajectory has large defects.

Model residual versus integration error

A zero defect for a coarse discrete scheme does not prove accuracy relative to the continuous model. Dynamic feasibility must be paired with time-step or mesh-convergence studies. Discontinuities, saturation transitions, impacts, and fast electrical modes may require local refinement or separate phases.

import numpy as np

def positioner_rhs(x, voltage, inertia=0.08, damping=0.12,
                   resistance=2.0, inductance=0.08,
                   torque_constant=0.25, back_emf=0.25):
    q, omega, current = x
    return np.array([
        omega,
        (torque_constant * current - damping * omega) / inertia,
        (voltage - resistance * current - back_emf * omega) / inductance,
    ])

dt = 0.01
x0 = np.array([0.0, 0.0, 0.0])
u0 = 4.0
x1_consistent = x0 + dt * positioner_rhs(x0, u0)
x1_inconsistent = x1_consistent + np.array([0.0, 0.05, 0.0])

def euler_defect(x0, x1, u0):
    return x1 - x0 - dt * positioner_rhs(x0, u0)

print("consistent defect:", euler_defect(x0, x1_consistent, u0))
print("inconsistent defect:", euler_defect(x0, x1_inconsistent, u0))

Activity 4.5: interpret a defect