Struct LqRegulator 

pub struct LqRegulator<const N: usize, const M: usize> {
    pub a: Vec<f64>,
    pub b: Vec<f64>,
    pub c: Vec<f64>,
    pub q: Vec<f64>,
    pub q_terminal: Vec<f64>,
    pub r: Vec<f64>,
    pub horizon: f64,
    /* private fields */
}

Finite-horizon linear-quadratic regulator with constant coefficients.

Matrices are stored row-major as Vec<f64>. The control dimension M is part of the type so the Control associated type can be a fixed-length array.

See the module-level documentation for the formulation and exact solution.

Fields

a: Vec<f64>

State-drift coupling matrix A (row-major, n x n).

b: Vec<f64>

Control-input matrix B (row-major, n x m).

c: Vec<f64>

Diffusion matrix C (row-major, n x n).

q: Vec<f64>

Running state cost Q (row-major, n x n).

q_terminal: Vec<f64>

Terminal state cost Q_T (row-major, n x n).

r: Vec<f64>

Running control cost R (row-major, m x m).

horizon: f64

Remaining horizon T.

Implementations

impl<const N: usize, const M: usize> LqRegulator<N, M>

pub fn new( a: &[f64], b: &[f64], c: &[f64], q: &[f64], q_terminal: &[f64], r: &[f64], horizon: f64, ) -> Self

Creates an LQ regulator.

Matrices are row-major. All matrix sizes are asserted to match N and M.

Panics

Panics if horizon is not positive or any matrix has the wrong length.

Examples

use solver::models::lq_regulator::LqRegulator;
let lq = LqRegulator::<2, 1>::new(
    &[0.0, 1.0, 0.0, 0.0],
    &[0.0, 1.0],
    &[0.0, 0.0, 0.0, 0.0],
    &[1.0, 0.0, 0.0, 0.0],
    &[1.0, 0.0, 0.0, 0.0],
    &[1.0],
    1.0,
);
assert_eq!(lq.horizon, 1.0);

pub fn riccati_solution(&self) -> (Vec<f64>, f64)

Solves the backward Riccati ODE with an implicit Euler scheme.

Returns (p, q) where p is the flat row-major solution of P(t) at t = 0 and q is the scalar noise-correction term q(0).

The number of time steps is fixed by Self::riccati_steps.

Examples

use solver::models::lq_regulator::LqRegulator;
let lq = LqRegulator::<1, 1>::new(
    &[0.0], &[1.0], &[0.0], &[1.0], &[1.0], &[1.0], 1.0,
);
let (p, q) = lq.riccati_solution();
assert!(p[0] > 0.0);
assert!(q >= 0.0);

pub fn riccati_solution_with_steps(&self, steps: usize) -> (Vec<f64>, f64)

Solves the backward Riccati ODE for the first steps fixed-size backward steps from the terminal condition.

Each step has size horizon / riccati_steps, so steps selects the remaining horizon steps * horizon / riccati_steps at which P and q are evaluated. steps = 0 returns the terminal condition (Q_T, 0); steps = riccati_steps() returns the solution at the start of the horizon.

pub fn riccati_steps(&self) -> usize

Number of Riccati backward steps.

The Riccati ODE is independent of the spatial grid and, for constant coefficients, its solution is essentially independent of the step size once dt is small. A fixed modest step count keeps the reference cheap and deterministic.

pub fn feedback_gain(&self, t: f64) -> Vec<f64>

Optimal feedback gain K = -R^{-1} B' P(t) at forward time t.

Returns a row-major m x n matrix K so that u = K x.

pub fn exact_value(&self, state: &[f64; N], tau: f64) -> f64

Exact value V(t, x) = -x' P(t) x - q(t) for remaining horizon tau.

pub fn exact_control(&self, state: &[f64; N], tau: f64) -> [f64; M]

Exact control u* = -R^{-1} B' P(t) x for remaining horizon tau.

Trait Implementations

impl<const N: usize, const M: usize> Clone for LqRegulator<N, M>

fn clone(&self) -> LqRegulator<N, M>

Returns a duplicate of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<const N: usize, const M: usize> ControlProblem<N> for LqRegulator<N, M>

type Control = [f64; M]

The action type. For Merton this is a scalar portfolio fraction; for market making it is a pair of bid/ask intensities.

fn optimize( &self, t: f64, state: &[f64; N], _derivs: &StateDerivatives<N>, ) -> Self::Control

Returns the control that maximizes the driver at (t, state).

fn running_reward( &self, _t: f64, state: &[f64; N], control: &Self::Control, ) -> f64

Returns the running reward f(t,x,u).

fn generator( &self, _t: f64, state: &[f64; N], control: &Self::Control, derivs: &StateDerivatives<N>, ) -> f64

Returns the infinitesimal generator L^u V for the given control and derivative bundle.

fn terminal(&self, state: &[f64; N]) -> f64

Terminal value g(x) at the horizon.

fn next_step( &self, t: f64, state: &[f64; N], dt: f64, noise: &[f64; N], ) -> [f64; N]

Advances the state one step under the optimal control at forward time t.

fn is_diffusion_dimension(&self, _dim: usize) -> bool

Whether a dimension is driven by Brownian diffusion.

fn driver( &self, t: f64, state: &[f64; N], control: &Self::Control, derivs: &StateDerivatives<N>, ) -> f64

Returns the full HJB driver f(t,x,u) + L^u V.

fn bsde_driver( &self, t: f64, state: &[f64; N], control: &Self::Control, _derivs: &StateDerivatives<N>, _dt: f64, ) -> f64

Returns the backward driver consumed by the BSDE regression solver.

fn apply_constraint(&self, _state: &[f64; N], value: f64) -> f64

Optional pointwise constraint on the value, e.g. the early-exercise obstacle V >= payoff for an American option.

fn discount_rate(&self, _state: &[f64; N]) -> f64

Discount rate r(t,x). Defaults to zero.

fn constant_discount_rate(&self) -> Option<f64>

Constant discount rate hint. Defaults to None (state dependent).

fn next_step_controlled( &self, t: f64, state: &[f64; N], _control: &Self::Control, dt: f64, noise: &[f64; N], ) -> [f64; N]

Advances the state one step under an explicitly supplied control.

fn is_reduced_value(&self) -> bool

Whether the value this problem solves is a reduced value.

fn gradient_step(&self, _dim: usize) -> f64

Physical finite-difference step for a dimension, used by mesh-free gradient stencils.

impl<const N: usize, const M: usize> Debug for LqRegulator<N, M>

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.