Incremental Model-based Global Dual Heuristic Programming (IMGDHP)¶

IMGDHP is an incremental model-based variant of Global Dual Heuristic Programming from the Adaptive Critic Designs (ACD) family. It is designed for online adaptive control of nonlinear systems under partial observability. The agent combines recursive least squares (RLS) system identification with a dual-head critic that estimates both the cost-to-go \(J\) and the costate vector \(\lambda = \partial J / \partial y\), enabling richer gradient information for the actor. See also the nonlinear F-16 model: NonlinearLongitudinalF16.

Key ideas¶

Incremental model: online identification of local linearization \(\Delta y_{t+1} = A \Delta y_t + B \Delta u_t\) via RLS — lightweight, interpretable, and does not require a neural network for system ID
GDHP dual critic: the critic outputs both \(J(o)\) (scalar cost-to-go) and \(\lambda(o)\) (costate vector), providing a richer gradient signal to the actor compared to standard HDP/DHP
Model-predictive actor update: the actor gradient flows through the identified model matrices \(A\), \(B\), enabling one-step lookahead optimization
Partial observability: the augmented observation \(o = [y; r; e]\) allows the agent to operate when the environment observation is not the full state

Aspect	HDP	IHDP	IMGDHP
System ID	Fixed/known model	Online NN	Online RLS (incremental linear)
Critic output	\(J(o)\) only	\(J(o)\) only	\(J(o)\) + \(\lambda(o)\) (dual)
Actor update	Direct gradient	Model-based	Model-predictive via \(A\), \(B\)
Partial observability	No	Limited	Core design feature
Framework	NumPy	NumPy	PyTorch

IMGDHP components¶

Component	Role	Implementation
IncrementalModelRLS	Online identification of \(A\), \(B\) matrices via RLS	`tensoraerospace.agent.im_gdhp.IncrementalModelRLS`
GDHPActor	Deterministic policy network \(u = u_{\max} \tanh(\pi_\theta(o))\)	`tensoraerospace.agent.im_gdhp.GDHPActor`
GDHPCritic	Dual-head critic: shared backbone with \(J\)-head and \(\lambda\)-head	`tensoraerospace.agent.im_gdhp.GDHPCritic`
IMGDHPAgent	Orchestrates all components, training loop, predict/learn interface	`tensoraerospace.agent.im_gdhp.IMGDHPAgent`

Algorithm¶

At each time step \(t\), given observation \(y_t\) and reference \(r_t\):

Augment observation: \(o_t = [y_t;\; r_t;\; e_t]\), where \(e_t = y_t[\text{tracking}] - r_t\)
Actor produces action: \(u_t = \pi_\theta(o_t)\)
Execute \(u_t\) in the environment, observe \(y_{t+1}\)
Compute one-step cost: \(c_t = e_t^\top Q e_t + \rho \| u_t - u_{t-1} \|^2\)
RLS update (if \(t \geq 2\)): update incremental model using \((y_{t-2}, y_{t-1}, y_t, u_{t-2}, u_{t-1})\) to obtain \(A_t\), \(B_t\)
Critic update (GDHP dual loss):

\[ L = \underbrace{\left( J(o_t) - (c_t + \gamma J(o_{t+1})) \right)^2}_{L_J} + \beta \underbrace{\left\| \lambda(o_t) - \left( \frac{\partial c_t}{\partial y} + \gamma A_t^\top \lambda(o_{t+1}) \right) \right\|^2}_{L_\lambda} \]

Actor update (model-predictive):

\[ \min_\theta \; c_t + \gamma \, J\!\left(\hat{o}_{t+1}\right), \quad \text{gradient flows through } B_t \]

Quick start¶

import numpy as np
import gymnasium as gym
from tensoraerospace.agent.im_gdhp import IMGDHPAgent, IMGDHPConfig
from tensoraerospace.utils import generate_time_period
from tensoraerospace.signals.standard import sinusoid

dt = 0.01
tp = generate_time_period(tn=20, dt=dt)
number_time_steps = len(tp)
reference_signal = sinusoid(
    degree=3, tp=tp, frequency=0.1, output_rad=True
).reshape(1, -1)

env = gym.make(
    "NonlinearLongitudinalF16-v0",
    number_time_steps=number_time_steps,
    initial_state=np.array([0.0, 0.0]),
    reference_signal=reference_signal,
    dt=dt,
)

config = IMGDHPConfig(
    gamma=0.95,
    actor_hidden=(32, 32),
    critic_hidden=(64, 64),
    actor_lr=1e-3,
    critic_lr=5e-3,
    track_Q=(1.0,),
    warmup_steps=5,
    forgetting=0.9995,
    u_max=25.0,
)

agent = IMGDHPAgent(
    n_obs=2,
    n_action=1,
    reference_size=1,
    tracking_indices=[0],
    config=config,
)

obs, info = env.reset()
for t in range(number_time_steps - 1):
    action = agent.predict(obs, reference_signal, t)
    obs_next, reward, terminated, truncated, info = env.step(action)
    metrics = agent.learn(obs_next, reference_signal, t)
    obs = obs_next
    if terminated or truncated:
        break

Tip

tracking_indices must align with the observation indices that correspond to the tracked reference signal. For example, if the observation is [alpha, wz] and you track alpha, use tracking_indices=[0].

Hyperparameters¶

General¶

Parameter	Default	Description
`gamma`	0.95	Discount factor
`warmup_steps`	5	Steps with frozen actor/critic (exploration only)
`critic_only_steps`	0	Additional steps with frozen actor after warmup
`seed`	None	RNG seed for reproducibility
`device`	`"cpu"`	PyTorch device

Actor¶

Parameter	Default	Description
`actor_hidden`	(32, 32)	Hidden layer sizes
`actor_lr`	1e-3	Learning rate
`u_max`	25.0	Per-channel control bound
`exploration_noise_std`	0.0	Gaussian exploration noise during training

Critic¶

Parameter	Default	Description
`critic_hidden`	(64, 64)	Backbone hidden layer sizes
`critic_lr`	5e-3	Learning rate
`beta_lambda`	1.0	Weight of \(\lambda\)-loss in GDHP dual objective
`critic_updates_per_step`	1	Gradient steps per environment transition
`target_update_tau`	0.0	Polyak coefficient for target critic (0 = no target network)
`critic_weight_decay`	0.0	L2 regularization
`max_grad_norm`	5.0	Gradient clipping threshold

Cost function¶

Parameter	Default	Description
`track_Q`	(1.0,)	Diagonal weights of tracking cost \(e^\top Q e\)
`action_rate_penalty`	1e-3	Coefficient \(\rho\) for \(\\| \Delta u \\|^2\) penalty

Incremental model (RLS)¶

Parameter	Default	Description
`forgetting`	0.9995	RLS forgetting factor \(\in (0, 1]\)
`cov_init`	1e2	Initial covariance matrix scale

Observation¶

Parameter	Default	Description
`obs_scale`	None	Per-component observation scaling factors

Supported environments¶

NonlinearLongitudinalF16-v0
LinearLongitudinalF16-v0

API reference¶

`IMGDHPAgent(n_obs, n_action, reference_size=1, tracking_indices=None, config=None)` ¶

Online IM-GDHP control agent with partial observability support.

The agent's public training interface mirrors the online style of :class:tensoraerospace.agent.ihdp.IHDPAgent: at each environment step the caller invokes :meth:predict with the latest observation and reference, gets back the next control command, executes it in the environment, then calls :meth:learn with the new observation to perform one RLS + critic + actor update. A convenience :meth:train loop wraps this for episodic training against a Gymnasium environment.

Parameters:

Name	Type	Description	Default
`n_obs`	`int`	Length of the environment observation vector `y`.	required
`n_action`	`int`	Number of control channels.	required
`reference_size`	`int`	Length of the reference vector at each time step (typically 1 per tracked channel).	`1`
`tracking_indices`	`Sequence[int] \| None`	Indices into `y` of the states that should track the reference. Used to build the scalar tracking error that drives the reward. Defaults to `[0]`.	`None`
`config`	`IMGDHPConfig \| None`	Optional :class:`IMGDHPConfig` instance. Fields not set fall back to the class defaults.	`None`

`reset()` ¶

Reset the per-episode rolling history.

Does not reset the learned weights, the incremental model covariance or the optimiser state — so learning progresses across episodes as expected.

`predict(obs, reference_signal, time_step, *, deterministic=False)` ¶

Compute the control action for a single time step.

Parameters:

Name	Type	Description	Default
`obs`	`ndarray`	Current observation `y_t` of length `n_obs`.	required
`reference_signal`	`ndarray`	Reference trajectory, shape `(reference_size, T)` or `(T,)`.	required
`time_step`	`int`	Current time index used to select `reference_signal[:, time_step]`.	required
`deterministic`	`bool`	If True, no exploration noise is added.	`False`

Returns:

Type	Description
`ndarray`	A NumPy array of length `n_action` with the commanded
`ndarray`	control.

`learn(next_obs, reference_signal, time_step)` ¶

Perform one online RLS + critic + actor update.

Must be called after :meth:predict and the corresponding environment step. next_obs is y_{t+1}, and time_step is the index of the step that has just been executed (i.e. the same time_step that was passed to :meth:predict).

Returns:

Type	Description
`dict[str, float]`	Dict with latest scalar training metrics (`critic_loss`,
`dict[str, float]`	`actor_loss`, `rls_pred_error_norm`).

`train(env, num_episodes=1, *, max_steps=None, verbose=False)` ¶

Run episodic training against a Gymnasium environment.

The env is expected to expose a reference_signal attribute (shape (reference_size, T)) — matching the convention used by the existing tensoraerospace.envs.f16 envs.

Parameters:

Name	Type	Description	Default
`env`	`Any`	Gymnasium-like environment to train on.	required
`num_episodes`	`int`	Number of episodes to run.	`1`
`max_steps`	`int \| None`	Optional cap on steps per episode.	`None`
`verbose`	`bool`	If True, print per-episode summaries.	`False`

Returns:

Type	Description
`dict[str, list[float]]`	The accumulated training history (same as `self.history`).

`get_param_env()` ¶

Build a JSON-serialisable config for :meth:save.

The agent has no bound environment (observations are fed in by the caller), so only the constructor signature and config dataclass are persisted.

`save(path=None, *, save_gradients=False)` ¶

Write the agent to a directory.

Files produced

config.json — constructor kwargs + serialised :class:IMGDHPConfig.
actor.pth / critic.pth / target_critic.pth — PyTorch state dicts for the three networks.
incremental_model.npz — RLS theta and P matrices plus scalar hyper-parameters.
actor_optim.pth / critic_optim.pth — optimiser state dicts (only when save_gradients=True).

Parameters:

Name	Type	Description	Default
`path`	`Union[str, Path, None]`	Base directory. If `None`, uses CWD.	`None`
`save_gradients`	`bool`	Persist optimiser states so training can resume bitwise-identically from the checkpoint.	`False`

Returns:

Type	Description
`str`	Absolute path to the created run directory.

`from_pretrained(repo_name, access_token=None, version=None, load_gradients=False)` `classmethod` ¶

Load an agent from a local directory or Hugging Face Hub.

Parameters:

Name	Type	Description	Default
`repo_name`	`str`	Local folder path, or `namespace/repo_name` on the Hugging Face Hub.	required
`access_token`	`Optional[str]`	Hub access token for private repos.	`None`
`version`	`Optional[str]`	Hub revision / branch / tag.	`None`
`load_gradients`	`bool`	Also restore optimiser state dicts.	`False`

Returns:

Name	Type	Description
`IMGDHPAgent`	`'IMGDHPAgent'`	Reconstructed agent.

`publish_to_hub(repo_name, folder_path, access_token=None)` ¶

Upload a :meth:save directory to the Hugging Face Hub.

Parameters:

Name	Type	Description	Default
`repo_name`	`str`	Target repository id, e.g. `"me/my-imgdhp"`.	required
`folder_path`	`Union[str, Path]`	Local folder produced by :meth:`save`.	required
`access_token`	`Optional[str]`	Hub access token.	`None`

`IMGDHPConfig(gamma=0.95, actor_hidden=(32, 32), critic_hidden=(64, 64), actor_lr=0.001, critic_lr=0.005, beta_lambda=1.0, track_Q=(1.0,), action_rate_penalty=0.001, forgetting=0.9995, cov_init=100.0, warmup_steps=5, critic_only_steps=0, critic_updates_per_step=1, target_update_tau=0.0, critic_weight_decay=0.0, obs_scale=None, max_grad_norm=5.0, exploration_noise_std=0.0, u_max=25.0, device='cpu', seed=None, history=dict())` `dataclass` ¶

Hyper-parameters for :class:IMGDHPAgent.

Parameters:

Name	Type	Description	Default
`gamma`	`float`	Discount factor γ.	`0.95`
`actor_hidden`	`Sequence[int]`	Hidden layer sizes of the actor MLP.	`(32, 32)`
`critic_hidden`	`Sequence[int]`	Hidden layer sizes of the critic backbone.	`(64, 64)`
`actor_lr`	`float`	Actor learning rate.	`0.001`
`critic_lr`	`float`	Critic learning rate.	`0.005`
`beta_lambda`	`float`	Weight of the λ-loss in the GDHP critic objective.	`1.0`
`track_Q`	`Sequence[float]`	Diagonal weights of the quadratic tracking cost over the tracked output channels. Length must equal `n_track`.	`(1.0,)`
`action_rate_penalty`	`float`	ρ coefficient penalising ‖Δu‖² in the one-step cost.	`0.001`
`forgetting`	`float`	RLS forgetting factor for the incremental model.	`0.9995`
`cov_init`	`float`	Initial scale of the RLS covariance matrix.	`100.0`
`warmup_steps`	`int`	Number of initial steps during which only the incremental model is updated (actor and critic are held fixed) so that `A` and `B` stabilise before being used for policy improvement.	`5`
`critic_only_steps`	`int`	Number of additional steps beyond `warmup_steps` during which the critic is updated but the actor is held fixed. This gives the value function time to settle before the policy gradient starts acting on it, which empirically dampens DHP oscillations.	`0`
`critic_updates_per_step`	`int`	Number of gradient steps to take on the critic for every environment transition. Multiple critic steps per env step accelerate Bellman convergence without waiting for more samples.	`1`
`target_update_tau`	`float`	Polyak coefficient for the soft target critic update `θ_target ← (1-τ)·θ_target + τ·θ`. When `0.0` the target network is disabled and the online critic is used for Bellman bootstrapping (legacy behaviour). Typical values `1e-3` to `1e-2`.	`0.0`
`critic_weight_decay`	`float`	L2 regularisation coefficient passed to the critic's `Adam` optimiser. Prevents the λ-head from drifting during long online runs.	`0.0`
`obs_scale`	`Sequence[float] \| None`	Optional per-component multiplier applied to the raw observation and the reference before feeding them to the actor and critic. Useful to compensate for very small (rad-scale) physical units. Length `n_obs` or `None`.	`None`
`max_grad_norm`	`float`	Gradient clipping applied to both actor and critic optimisers.	`5.0`
`exploration_noise_std`	`float`	Std of zero-mean Gaussian noise added to the actor output during training. Set to `0` to disable.	`0.0`
`u_max`	`float`	Per-channel absolute bound on the control output of the actor, in the units expected by the environment (e.g. deg for the F-16 envs in tensoraerospace).	`25.0`
`device`	`str`	Torch device for the networks.	`'cpu'`
`seed`	`int \| None`	Optional seed for torch / numpy RNGs.	`None`

`IncrementalModelRLS(n_y, n_u, forgetting=0.999, cov_init=100.0, theta_init_scale=0.001, seed=None)` ¶

Recursive least squares identifier for the incremental model.

Maintains a parameter matrix theta of shape (n_y + n_u, n_y) such that the first n_y rows correspond to Aᵀ and the last n_u rows correspond to Bᵀ. The identification equation is::

Δy_{t+1}ᵀ ≈ φ_tᵀ · theta

with regressor φ_t = [Δy_t; Δu_t] of length n_y + n_u.

Parameters:

Name	Type	Description	Default
`n_y`	`int`	Dimension of the observed/tracked output vector y.	required
`n_u`	`int`	Dimension of the control input vector u.	required
`forgetting`	`float`	RLS forgetting factor α ∈ (0, 1]. Values <1 give more weight to recent samples; `1.0` recovers ordinary least squares. Typical flight-control values: 0.995–0.9999.	`0.999`
`cov_init`	`float`	Initial scale of the covariance matrix P₀ = cov_init · I. Larger values encourage faster initial learning at the cost of robustness to noisy early samples.	`100.0`
`theta_init_scale`	`float`	Standard deviation of the zero-mean Gaussian used to randomly initialise `theta`. A small non-zero value helps break the symmetry when later matrix operations are used (e.g. `Bᵀ B`).	`0.001`
`seed`	`int \| None`	Optional random seed for `theta` initialisation.	`None`

`A` `property` ¶

Return the identified A matrix (shape (n_y, n_y)).

`B` `property` ¶

Return the identified B matrix (shape (n_y, n_u)).

`reset()` ¶

Reset the parameter buffers (keeps theta and P).

`reset_covariance()` ¶

Reset P to its initial large-variance state (full re-learn).

`update(y_prev, y_curr, y_next, u_prev, u_curr)` ¶

Perform a single RLS step from a sliding window of length 3.

Uses the tuple (y_{t-1}, y_t, y_{t+1}, u_{t-1}, u_t) to form the regressor φ = [y_t − y_{t-1}; u_t − u_{t-1}] and the target δ = y_{t+1} − y_t.

Parameters:

Name	Type	Description	Default
`y_prev`	`ndarray`	`y_{t-1}` — observation two steps back.	required
`y_curr`	`ndarray`	`y_t` — observation one step back.	required
`y_next`	`ndarray`	`y_{t+1}` — most recent observation.	required
`u_prev`	`ndarray`	`u_{t-1}` — control applied at step `t-1`.	required
`u_curr`	`ndarray`	`u_t` — control applied at step `t`.	required

Returns:

Type	Description
`ndarray`	The prediction error `ε = δ − θᵀ φ` before the update.

`predict_next(y_curr, y_prev, u_curr, u_prev)` ¶

Predict y_{t+1} from the current model estimates.

Uses the incremental form y_{t+1} ≈ y_t + A · (y_t − y_{t-1}) +B · (u_t − u_{t-1}).

`GDHPActor(in_features, n_u, hidden_sizes=(32, 32), u_max=1.0, activation=nn.Tanh)` ¶

Bases: Module

Deterministic policy network for IM-GDHP.

Maps an augmented observation [y; y_ref; e_track] to a bounded control command u ∈ [-u_max, u_max]. Bounding is enforced by a tanh output head scaled by u_max — this matches the "symmetrical sigmoid activation in the output layer of the actor" trick from the Sun/van Kampen incremental ADP papers.

Parameters:

Name	Type	Description	Default
`in_features`	`int`	Size of the augmented observation vector.	required
`n_u`	`int`	Number of control channels.	required
`hidden_sizes`	`Sequence[int]`	Sizes of the hidden layers.	`(32, 32)`
`u_max`	`float`	Per-channel absolute bound on the control command.	`1.0`
`activation`	`type[Module]`	Hidden-layer activation factory.	`Tanh`

`forward(obs)` ¶

Compute u = u_max · tanh(head(backbone(obs))).

`GDHPCritic(in_features, n_y, hidden_sizes=(32, 32), activation=nn.Tanh)` ¶

Bases: Module

Dual-head critic for Global Dual Heuristic Programming.

Outputs both the scalar cost-to-go J(o) and its vector derivative λ(o) = ∂J/∂y (same dimensionality as the observed state y). The two heads share a common backbone so that the J-regression and the λ-regression reinforce each other during training.

Parameters:

Name	Type	Description	Default
`in_features`	`int`	Size of the augmented observation vector `o` (typically `n_y + n_ref + n_y` when the tracking error is concatenated).	required
`n_y`	`int`	Size of the observed state vector `y`. Defines the number of outputs of the λ head.	required
`hidden_sizes`	`Sequence[int]`	Sizes of the shared hidden layers.	`(32, 32)`
`activation`	`type[Module]`	Hidden-layer activation factory.	`Tanh`

`forward(obs)` ¶

Return (J, λ) for a batch of augmented observations.

Sources¶

Sun, Z. & van Kampen, E.-J. (2021). Intelligent adaptive optimal control using incremental model-based global dual heuristic programming subject to partial observability. Applied Soft Computing, 103, 107153.
Zhou, Y., van Kampen, E.-J., & Chu, Q. P. (2020). Incremental model based online dual heuristic programming for nonlinear adaptive control. Control Engineering Practice, 95, 104242.

Incremental Model-based Global Dual Heuristic Programming (IMGDHP)¶

Key ideas¶

Key differences from related agents¶

IMGDHP components¶

Algorithm¶

Quick start¶

Hyperparameters¶

General¶

Actor¶

Critic¶

Cost function¶

Incremental model (RLS)¶

Observation¶

Supported environments¶

API reference¶

IMGDHPAgent(n_obs, n_action, reference_size=1, tracking_indices=None, config=None) ¶

reset() ¶

predict(obs, reference_signal, time_step, *, deterministic=False) ¶

learn(next_obs, reference_signal, time_step) ¶

train(env, num_episodes=1, *, max_steps=None, verbose=False) ¶

get_param_env() ¶

save(path=None, *, save_gradients=False) ¶

from_pretrained(repo_name, access_token=None, version=None, load_gradients=False) classmethod ¶

publish_to_hub(repo_name, folder_path, access_token=None) ¶

IncrementalModelRLS(n_y, n_u, forgetting=0.999, cov_init=100.0, theta_init_scale=0.001, seed=None) ¶

A property ¶

B property ¶

reset() ¶

reset_covariance() ¶

update(y_prev, y_curr, y_next, u_prev, u_curr) ¶

predict_next(y_curr, y_prev, u_curr, u_prev) ¶

GDHPActor(in_features, n_u, hidden_sizes=(32, 32), u_max=1.0, activation=nn.Tanh) ¶

forward(obs) ¶

GDHPCritic(in_features, n_y, hidden_sizes=(32, 32), activation=nn.Tanh) ¶

forward(obs) ¶

Sources¶

`IMGDHPAgent(n_obs, n_action, reference_size=1, tracking_indices=None, config=None)` ¶

`reset()` ¶

`predict(obs, reference_signal, time_step, *, deterministic=False)` ¶

`learn(next_obs, reference_signal, time_step)` ¶

`train(env, num_episodes=1, *, max_steps=None, verbose=False)` ¶

`get_param_env()` ¶

`save(path=None, *, save_gradients=False)` ¶

`from_pretrained(repo_name, access_token=None, version=None, load_gradients=False)` `classmethod` ¶

`publish_to_hub(repo_name, folder_path, access_token=None)` ¶

`IncrementalModelRLS(n_y, n_u, forgetting=0.999, cov_init=100.0, theta_init_scale=0.001, seed=None)` ¶

`A` `property` ¶

`B` `property` ¶

`reset()` ¶

`reset_covariance()` ¶

`update(y_prev, y_curr, y_next, u_prev, u_curr)` ¶

`predict_next(y_curr, y_prev, u_curr, u_prev)` ¶

`GDHPActor(in_features, n_u, hidden_sizes=(32, 32), u_max=1.0, activation=nn.Tanh)` ¶

`forward(obs)` ¶

`GDHPCritic(in_features, n_y, hidden_sizes=(32, 32), activation=nn.Tanh)` ¶

`forward(obs)` ¶