optimalcontrol

package
v0.8.0 Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: Jul 19, 2026 License: MIT Imports: 5 Imported by: 0

Documentation

Overview

Package optimalcontrol implements optimal control and dynamic programming algorithms using only the Go standard library.

The package is self-contained: it ships its own small dense linear-algebra layer (the Matrix type with LU, Cholesky, least-squares, matrix exponential, Faddeev–LeVerrier characteristic polynomials and Durand–Kerner eigenvalues, symmetric Jacobi eigensolves, and Lyapunov/Sylvester solvers) and builds the control-theoretic algorithms on top of it.

The control content spans the classical pillars of the field:

  • Linear-quadratic regulators in continuous and discrete time, obtained by solving the algebraic Riccati equations. The continuous CARE is solved by the matrix-sign-function method applied to the Hamiltonian matrix and by Kleinman's Newton iteration; the discrete DARE is solved by the Riccati recursion. See SolveCARE, SolveDARE, LQRContinuous and LQRDiscrete.

  • Finite-horizon control via the backward Riccati recursion (discrete) and backward integration of the matrix Riccati differential equation (continuous). See FiniteHorizonLQRDiscrete and FiniteHorizonLQRContinuous.

  • Pontryagin's minimum principle expressed as a Hamiltonian two-point boundary-value problem, solved exactly for linear-quadratic problems and by indirect single shooting for nonlinear ones. See PontryaginLQ and IndirectShooting.

  • Hamilton–Jacobi–Bellman value iteration on state grids via a semi-Lagrangian scheme. See HJBGrid1D.

  • Dynamic programming for finite Markov decision processes: value iteration, Gauss–Seidel value iteration, exact and iterative policy evaluation, Howard's policy iteration and modified policy iteration. See the MDP type.

  • Kalman filtering and the linear-quadratic-Gaussian dual of the LQR, including steady-state continuous and discrete filters and a recursive KalmanFilter. See KalmanContinuous, KalmanDiscrete, LQGContinuous and LQGDiscrete.

Structural analysis (controllability, observability, stabilizability and detectability via the Popov–Belevitch–Hautus tests) and Lyapunov Gramians round out the toolkit.

All algorithms are deterministic; where randomness is useful the caller supplies a seed and uses math/rand. Nothing outside the standard library is imported.

Index

Examples

Constants

This section is empty.

Variables

View Source 
var ErrDim = errors.New("optimalcontrol: incompatible dimensions")

ErrDim is returned when matrix or vector dimensions are incompatible.

View Source 
var ErrNotConverged = errors.New("optimalcontrol: iteration did not converge")

ErrNotConverged is returned by iterative solvers that fail to converge within the allotted number of iterations.

View Source 
var ErrSingular = errors.New("optimalcontrol: matrix is singular")

ErrSingular is returned when a matrix that must be invertible is singular.

Functions

func CharPoly 

func CharPoly(a *Matrix) []float64

CharPoly returns the coefficients of the characteristic polynomial det(λI − A) of a square matrix, in ascending powers of λ (index i holds the coefficient of λ^i). It uses the Faddeev–LeVerrier algorithm.

func ClosedLoopEigenvaluesContinuous 

func ClosedLoopEigenvaluesContinuous(a, b, k *Matrix) []complex128

ClosedLoopEigenvaluesContinuous returns the eigenvalues of A − B K, the continuous closed-loop dynamics under state feedback K.

func ClosedLoopEigenvaluesDiscrete 

func ClosedLoopEigenvaluesDiscrete(a, b, k *Matrix) []complex128

ClosedLoopEigenvaluesDiscrete returns the eigenvalues of A − B K, the discrete closed-loop dynamics under state feedback K.

func ControllabilityRank 

func ControllabilityRank(a, b *Matrix) int

ControllabilityRank returns the rank of the controllability matrix of (A, B).

func CostateDynamicsLQ 

func CostateDynamicsLQ(a, q *Matrix, x, p []float64) []float64

CostateDynamicsLQ returns the costate derivative p' = −Q x − Aᵀ p for the linear-quadratic problem.

func Det 

func Det(a *Matrix) (float64, error)

Det returns the determinant of a square matrix.

func Eigenvalues 

func Eigenvalues(a *Matrix) []complex128

Eigenvalues returns all (generally complex) eigenvalues of a square matrix by forming its characteristic polynomial and finding its roots with the Durand–Kerner method. Results are not returned in any particular order.

func IsControllable 

func IsControllable(a, b *Matrix) bool

IsControllable reports whether the pair (A, B) is controllable, i.e. the controllability matrix has full row rank n.

func IsDetectableContinuous 

func IsDetectableContinuous(a, c *Matrix) bool

IsDetectableContinuous reports whether the pair (A, C) is detectable in continuous time (dual of stabilizability).

func IsDetectableDiscrete 

func IsDetectableDiscrete(a, c *Matrix) bool

IsDetectableDiscrete reports whether the pair (A, C) is detectable in discrete time.

func IsObservable 

func IsObservable(a, c *Matrix) bool

IsObservable reports whether the pair (A, C) is observable, i.e. the observability matrix has full column rank n.

func IsPositiveDefinite 

func IsPositiveDefinite(a *Matrix) bool

IsPositiveDefinite reports whether a symmetric matrix is positive definite.

func IsPositiveSemiDefinite 

func IsPositiveSemiDefinite(a *Matrix, tol float64) bool

IsPositiveSemiDefinite reports whether a symmetric matrix has all eigenvalues >= -tol.

func IsStabilizableContinuous 

func IsStabilizableContinuous(a, b *Matrix) bool

IsStabilizableContinuous reports whether the pair (A, B) is stabilizable in continuous time: every eigenvalue of A with non-negative real part is controllable (Popov–Belevitch–Hautus test).

func IsStabilizableDiscrete 

func IsStabilizableDiscrete(a, b *Matrix) bool

IsStabilizableDiscrete reports whether the pair (A, B) is stabilizable in discrete time: every eigenvalue of A with modulus >= 1 is controllable.

func IsStableContinuous 

func IsStableContinuous(a *Matrix, tol float64) bool

IsStableContinuous reports whether every eigenvalue of A has real part below -tol (Hurwitz stability of x' = A x).

func IsStableDiscrete 

func IsStableDiscrete(a *Matrix, tol float64) bool

IsStableDiscrete reports whether every eigenvalue of A has modulus below 1-tol (Schur stability of x_{k+1} = A x_k).

func LQRCostToGo 

func LQRCostToGo(p *Matrix, x []float64) float64

LQRCostToGo returns the optimal quadratic cost xᵀ P x of steering the state x to the origin under an infinite-horizon LQR with Riccati solution P.

func LeastSquares 

func LeastSquares(a *Matrix, b []float64) ([]float64, error)

LeastSquares solves the overdetermined system A x = b in the least-squares sense via the normal equations AᵀA x = Aᵀb. A must have full column rank.

func LinearInterp1D 

func LinearInterp1D(grid, values []float64, x float64) float64

LinearInterp1D linearly interpolates the tabulated function (grid, values) at x. The grid must be sorted ascending. Queries outside the grid are clamped to the endpoint values.

func ObservabilityRank 

func ObservabilityRank(a, c *Matrix) int

ObservabilityRank returns the rank of the observability matrix of (A, C).

func OptimalControlLQ 

func OptimalControlLQ(b, r *Matrix, p []float64) ([]float64, error)

OptimalControlLQ returns the stationarity-condition control u = −R⁻¹ Bᵀ p that minimizes the linear-quadratic Hamiltonian for a given costate p.

func PolyRootsDK 

func PolyRootsDK(coeffs []float64) []complex128

PolyRootsDK finds all complex roots of a real polynomial using the Durand–Kerner (Weierstrass) iteration. Coefficients are given in ascending power order: coeffs[i] multiplies x^i. The leading coefficient must be nonzero.

func QuadraticCostDiscrete 

func QuadraticCostDiscrete(xs [][]float64, us [][]float64, q, r, qf *Matrix) float64

QuadraticCostDiscrete evaluates the discrete quadratic cost Σ_k (xₖᵀ Q xₖ + uₖᵀ R uₖ) + x_Nᵀ Qf x_N over a state trajectory and control sequence. The control sequence has one fewer entry than the state trajectory.

func Rank 

func Rank(a *Matrix, tol float64) int

Rank returns the numerical rank of a matrix using Gaussian elimination with partial pivoting and the supplied tolerance for treating pivots as zero.

func SimulateDiscreteLQR 

func SimulateDiscreteLQR(a, b *Matrix, fh *FiniteHorizonDiscrete, x0 []float64) [][]float64

SimulateDiscreteLQR simulates the closed-loop discrete system under the time-varying gains produced by FiniteHorizonLQRDiscrete, starting from x0 and returning the state trajectory x[0..N].

func Solve 

func Solve(a *Matrix, b []float64) ([]float64, error)

Solve solves the linear system A x = b.

func SpectralAbscissa 

func SpectralAbscissa(a *Matrix) float64

SpectralAbscissa returns the maximum real part of the eigenvalues of a matrix.

func SpectralRadius 

func SpectralRadius(a *Matrix) float64

SpectralRadius returns the maximum modulus of the eigenvalues of a matrix.

func SymEigenvalues 

func SymEigenvalues(a *Matrix) ([]float64, error)

SymEigenvalues returns the eigenvalues of a symmetric matrix in ascending order.

func VecAdd 

func VecAdd(a, b []float64) []float64

VecAdd returns the element-wise sum a+b.

func VecAxpy 

func VecAxpy(s float64, b, a []float64) []float64

VecAxpy returns a + s·b (the classic "axpy" operation).

func VecCopy 

func VecCopy(a []float64) []float64

VecCopy returns a copy of a.

func VecDot 

func VecDot(a, b []float64) float64

VecDot returns the inner product aᵀb.

func VecMaxAbs 

func VecMaxAbs(a []float64) float64

VecMaxAbs returns the largest absolute component of a.

func VecNorm 

func VecNorm(a []float64) float64

VecNorm returns the Euclidean norm of a.

func VecScale 

func VecScale(a []float64, s float64) []float64

VecScale returns the vector a scaled by s.

func VecSub 

func VecSub(a, b []float64) []float64

VecSub returns the element-wise difference a−b.

Types

type FiniteHorizonContinuous 

type FiniteHorizonContinuous struct {
	// Times holds the sample times, ascending from 0 to T.
	Times []float64
	// P holds the Riccati matrix at each sample time.
	P []*Matrix
}

FiniteHorizonContinuous holds the sampled solution of the continuous finite-horizon Riccati differential equation on a time grid.

func FiniteHorizonLQRContinuous 

func FiniteHorizonLQRContinuous(a, b, q, r, qf *Matrix, tFinal float64, steps int) (*FiniteHorizonContinuous, error)

FiniteHorizonLQRContinuous integrates the matrix Riccati differential equation 

−dP/dt = Aᵀ P + P A − P B R⁻¹ Bᵀ P + Q,   P(T) = Qf,

backward from the terminal time T to 0 using classical fourth-order Runge–Kutta with the given number of steps. Samples are returned in ascending time order.

func (*FiniteHorizonContinuous) GainAt 

func (fh *FiniteHorizonContinuous) GainAt(b, r *Matrix, i int) (*Matrix, error)

GainAt returns the continuous LQR feedback gain K(t) = R⁻¹ Bᵀ P(t) at the grid index i of a FiniteHorizonContinuous solution.

type FiniteHorizonDiscrete 

type FiniteHorizonDiscrete struct {
	// P holds the cost-to-go matrices P[0..N]; P[N] is the terminal weight.
	P []*Matrix
	// K holds the feedback gains K[0..N-1] with u_k = −K[k] x_k.
	K []*Matrix
}

FiniteHorizonDiscrete holds the time-varying gains and cost matrices of a finite-horizon discrete LQR problem, indexed from stage 0 (initial) to N (terminal).

func FiniteHorizonLQRDiscrete 

func FiniteHorizonLQRDiscrete(a, b, q, r, qf *Matrix, n int) (*FiniteHorizonDiscrete, error)

FiniteHorizonLQRDiscrete solves the finite-horizon discrete LQR problem for x_{k+1} = A x_k + B u_k minimizing xₙᵀ Qf xₙ + Σ_{k<N} (xₖᵀ Q xₖ + uₖᵀ R uₖ) by backward Riccati recursion. It returns the sequence of cost-to-go matrices and time-varying gains.

type HJBGrid1D 

type HJBGrid1D struct {
	// Grid is the ascending state grid.
	Grid []float64
	// Controls is the finite set of admissible control values searched at each
	// grid point.
	Controls []float64
	// Dynamics returns f(x, u), the state velocity.
	Dynamics func(x, u float64) float64
	// RunningCost returns L(x, u), the instantaneous cost.
	RunningCost func(x, u float64) float64
	// Dt is the time step of the semi-Lagrangian discretization.
	Dt float64
	// Rho is the (non-negative) discount rate.
	Rho float64
}

HJBGrid1D describes a one-dimensional infinite-horizon optimal-control problem discretized for solution by value iteration (a semi-Lagrangian scheme for the Hamilton–Jacobi–Bellman equation ρ V = min_u [ L(x,u) + V'(x) f(x,u) ]).

func (*HJBGrid1D) ControlAt 

func (h *HJBGrid1D) ControlAt(res *HJBResult, x float64) float64

ControlAt returns the greedy control at state x by interpolating the tabulated policy of a solved HJBResult onto the grid.

func (*HJBGrid1D) Solve 

func (h *HJBGrid1D) Solve(tol float64, maxIter int) *HJBResult

Solve runs value iteration for the infinite-horizon discounted problem until the max-norm change is below tol or maxIter sweeps elapse. The per-step discount is exp(−ρ Δt) and the successor value is obtained by linear interpolation of the current value function at x + Δt f(x, u).

func (*HJBGrid1D) SolveFiniteHorizon 

func (h *HJBGrid1D) SolveFiniteHorizon(stages int, terminal func(x float64) float64) *HJBResult

SolveFiniteHorizon solves the finite-horizon problem with terminal cost g(x) by backward value iteration over the given number of stages, returning the value function and control at the initial time (stage 0).

type HJBResult 

type HJBResult struct {
	// Value is the optimal cost-to-go at each grid point.
	Value []float64
	// Control is the optimal control at each grid point.
	Control []float64
	// Iterations is the number of value-iteration sweeps performed.
	Iterations int
	// Converged reports whether the sweep residual fell below the tolerance.
	Converged bool
}

HJBResult holds the value function and greedy control at each grid point.

type KalmanContinuousResult 

type KalmanContinuousResult struct {
	// P is the steady-state estimation-error covariance.
	P *Matrix
	// L is the Kalman gain L = P Cᵀ V⁻¹.
	L *Matrix
}

KalmanContinuousResult holds the steady-state error covariance and estimator gain of a continuous-time Kalman–Bucy filter.

func KalmanContinuous 

func KalmanContinuous(a, c, w, v *Matrix) (*KalmanContinuousResult, error)

KalmanContinuous designs the steady-state continuous-time Kalman–Bucy filter for x' = A x + w, y = C x + v with process covariance W and measurement covariance V. It solves the filter algebraic Riccati equation A P + P Aᵀ − P Cᵀ V⁻¹ C P + W = 0 (the dual of the LQR CARE) and returns the covariance and gain.

type KalmanDiscreteResult 

type KalmanDiscreteResult struct {
	// P is the steady-state a-priori (predicted) error covariance.
	P *Matrix
	// L is the steady-state Kalman update gain L = P Cᵀ (C P Cᵀ + V)⁻¹.
	L *Matrix
}

KalmanDiscreteResult holds the steady-state a-priori covariance and gain of a discrete-time Kalman filter.

func KalmanDiscrete 

func KalmanDiscrete(a, c, w, v *Matrix) (*KalmanDiscreteResult, error)

KalmanDiscrete designs the steady-state discrete-time Kalman filter for x_{k+1} = A x_k + w_k, y_k = C x_k + v_k with process covariance W and measurement covariance V. It solves the dual DARE P = A P Aᵀ − A P Cᵀ (C P Cᵀ + V)⁻¹ C P Aᵀ + W and returns the a-priori covariance and update gain.

type KalmanFilter 

type KalmanFilter struct {
	// A is the state-transition matrix.
	A *Matrix
	// B is the (optional) control-input matrix; may be nil for no input.
	B *Matrix
	// C is the measurement matrix.
	C *Matrix
	// Q is the process-noise covariance.
	Q *Matrix
	// R is the measurement-noise covariance.
	R *Matrix
	// X is the current state estimate.
	X []float64
	// P is the current estimate covariance.
	P *Matrix
}

KalmanFilter is a recursive discrete-time Kalman filter that tracks a state estimate and its covariance. Construct one with NewKalmanFilter and drive it with alternating Predict and Update calls.

func NewKalmanFilter 

func NewKalmanFilter(a, b, c, q, r *Matrix, x0 []float64, p0 *Matrix) *KalmanFilter

NewKalmanFilter constructs a Kalman filter with the given model matrices and initial estimate x0 with covariance p0. B may be nil when the system has no control input.

func (*KalmanFilter) Predict 

func (kf *KalmanFilter) Predict(u []float64)

Predict advances the estimate through the process model using control input u (which may be nil when B is nil): x̂ ← A x̂ + B u, P ← A P Aᵀ + Q.

func (*KalmanFilter) Update 

func (kf *KalmanFilter) Update(y []float64) error

Update corrects the prediction with measurement y using the Kalman gain: L = P Cᵀ (C P Cᵀ + R)⁻¹, x̂ ← x̂ + L (y − C x̂), P ← (I − L C) P. It returns an error if the innovation covariance is singular.

type LQGResult 

type LQGResult struct {
	// K is the LQR state-feedback gain (u = −K x̂).
	K *Matrix
	// L is the Kalman estimator gain.
	L *Matrix
	// P is the control Riccati solution.
	P *Matrix
	// Sigma is the estimation-error covariance.
	Sigma *Matrix
}

LQGResult bundles the regulator and estimator designs of a linear-quadratic Gaussian controller.

func LQGContinuous 

func LQGContinuous(a, b, c, q, r, w, v *Matrix) (*LQGResult, error)

LQGContinuous designs a continuous-time LQG controller by combining a continuous LQR (weights Q, R on the plant (A, B)) with a continuous Kalman filter (process covariance W, measurement covariance V on (A, C)), invoking the separation principle.

func LQGDiscrete 

func LQGDiscrete(a, b, c, q, r, w, v *Matrix) (*LQGResult, error)

LQGDiscrete designs a discrete-time LQG controller by combining a discrete LQR with a discrete Kalman filter via the separation principle.

type LQRResult 

type LQRResult struct {
	// K is the optimal feedback gain (u = −K x).
	K *Matrix
	// P is the stabilizing Riccati solution (the cost-to-go Hessian).
	P *Matrix
	// ClosedLoop is A − B K.
	ClosedLoop *Matrix
}

LQRResult bundles the gain matrix, Riccati solution and closed-loop matrix of a linear-quadratic regulator design.

func LQRContinuous 

func LQRContinuous(a, b, q, r *Matrix) (*LQRResult, error)

LQRContinuous designs a continuous-time infinite-horizon LQR for the system x' = A x + B u minimizing ∫ (xᵀ Q x + uᵀ R u) dt. The optimal control is u = −K x with K = R⁻¹ Bᵀ P and P the stabilizing CARE solution.

Example 
// Double integrator with unit weights: known gain [1, √3].
a := FromRows([][]float64{{0, 1}, {0, 0}})
b := FromRows([][]float64{{0}, {1}})
q := Eye(2)
r := Eye(1)
res, _ := LQRContinuous(a, b, q, r)
fmt.Printf("K = [%.4f %.4f]\n", res.K.At(0, 0), res.K.At(0, 1))

Output:
K = [1.0000 1.7321]

func LQRDiscrete 

func LQRDiscrete(a, b, q, r *Matrix) (*LQRResult, error)

LQRDiscrete designs a discrete-time infinite-horizon LQR for the system x_{k+1} = A x_k + B u_k minimizing Σ (xᵀ Q x + uᵀ R u). The optimal control is u = −K x with K = (R + Bᵀ P B)⁻¹ Bᵀ P A and P the stabilizing DARE solution.

type LU 

type LU struct {
	// contains filtered or unexported fields
}

LU holds the result of an LU factorization with partial pivoting: a combined lower/upper factor lu, the row permutation piv, and the sign of the permutation (+1 or -1) used for the determinant.

func Factor 

func Factor(a *Matrix) (*LU, error)

Factor computes the LU decomposition of a square matrix with partial pivoting. It returns ErrDim for non-square input.

func (*LU) Det 

func (f *LU) Det() float64

Det returns the determinant computed from the factorization.

func (*LU) Solve 

func (f *LU) Solve(b []float64) ([]float64, error)

Solve solves A x = b for x using the factorization.

func (*LU) SolveMatrix 

func (f *LU) SolveMatrix(b *Matrix) (*Matrix, error)

SolveMatrix solves A X = B, treating each column of B independently.

type LinearSystem 

type LinearSystem struct {
	// A is the n×n state matrix.
	A *Matrix
	// B is the n×m input matrix.
	B *Matrix
	// C is the p×n output matrix.
	C *Matrix
	// D is the p×m feedthrough matrix.
	D *Matrix
}

LinearSystem is a continuous- or discrete-time linear time-invariant model x' = A x + B u (continuous) or x_{k+1} = A x_k + B u_k (discrete), with output y = C x + D u.

func NewLinearSystem 

func NewLinearSystem(a, b, c, d *Matrix) *LinearSystem

NewLinearSystem builds a LinearSystem, filling C with the identity and D with zeros when they are nil.

func (*LinearSystem) Output 

func (s *LinearSystem) Output(x, u []float64) []float64

Output returns y = C x + D u for the system.

func (*LinearSystem) SimulateContinuous 

func (s *LinearSystem) SimulateContinuous(x0 []float64, control func(t float64) []float64, dt float64, steps int) [][]float64

SimulateContinuous integrates the continuous-time system with a fixed control u(t) held over each step of size dt using classical RK4, returning the state trajectory at the sample times 0, dt, …, steps·dt.

func (*LinearSystem) SimulateDiscrete 

func (s *LinearSystem) SimulateDiscrete(x0 []float64, us [][]float64) [][]float64

SimulateDiscrete simulates the discrete-time system for the given control sequence, returning the state trajectory x[0..len(us)].

func (*LinearSystem) StepDiscrete 

func (s *LinearSystem) StepDiscrete(x, u []float64) []float64

StepDiscrete advances the discrete-time state one step: x_{k+1} = A x + B u.

type MDP 

type MDP struct {
	// States is the number of states.
	States int
	// Actions is the number of actions.
	Actions int
	// Trans[a] is the States×States transition matrix under action a; each row
	// must sum to one.
	Trans []*Matrix
	// Reward is a States×Actions matrix of expected immediate rewards.
	Reward *Matrix
	// Gamma is the discount factor in [0, 1).
	Gamma float64
}

MDP is a finite, discounted Markov decision process with the standard reward-maximization convention. Transitions for each action are stored as row-stochastic matrices and rewards as expected immediate rewards per state-action pair.

func NewMDP 

func NewMDP(trans []*Matrix, reward *Matrix, gamma float64) (*MDP, error)

NewMDP constructs an MDP from per-action transition matrices, a reward matrix and a discount factor, validating the dimensions.

func (*MDP) BellmanBackup 

func (m *MDP) BellmanBackup(v []float64) (newV []float64, policy []int)

BellmanBackup applies one Bellman optimality update, returning the improved value function and the greedy policy that attains it.

func (*MDP) GaussSeidelValueIteration 

func (m *MDP) GaussSeidelValueIteration(tol float64, maxIter int) *ValueIterationResult

GaussSeidelValueIteration runs value iteration using in-place (Gauss–Seidel) updates, which typically converges in fewer sweeps than the Jacobi form.

func (*MDP) GreedyPolicy 

func (m *MDP) GreedyPolicy(v []float64) []int

GreedyPolicy returns the greedy policy with respect to a value function.

func (*MDP) ModifiedPolicyIteration 

func (m *MDP) ModifiedPolicyIteration(k, maxIter int, tol float64) (*PolicyIterationResult, error)

ModifiedPolicyIteration runs modified policy iteration, evaluating each policy with a fixed number k of Bellman expectation sweeps rather than an exact solve.

func (*MDP) PolicyEvaluationExact 

func (m *MDP) PolicyEvaluationExact(policy []int) ([]float64, error)

PolicyEvaluationExact computes the exact value function of a deterministic policy by solving the linear system (I − γ P_π) V = r_π.

func (*MDP) PolicyEvaluationIterative 

func (m *MDP) PolicyEvaluationIterative(policy []int, tol float64, maxIter int) []float64

PolicyEvaluationIterative computes the value function of a policy by iterative application of the Bellman expectation operator.

func (*MDP) PolicyIteration 

func (m *MDP) PolicyIteration(maxIter int) (*PolicyIterationResult, error)

PolicyIteration runs Howard's policy iteration with exact policy evaluation, returning the optimal policy and its value function.

func (*MDP) PolicyMatrices 

func (m *MDP) PolicyMatrices(policy []int) (*Matrix, []float64)

PolicyMatrices returns the transition matrix P_π and reward vector r_π induced by a deterministic policy.

func (*MDP) QValues 

func (m *MDP) QValues(v []float64) *Matrix

QValues returns the States×Actions matrix of action values Q(s, a) = R(s, a) + γ Σ_{s'} P_a(s, s') V(s') for a given value function.

func (*MDP) ValueIteration 

func (m *MDP) ValueIteration(tol float64, maxIter int) *ValueIterationResult

ValueIteration runs value iteration until the max-norm change between sweeps falls below tol or maxIter sweeps have elapsed.

Example 
stay := FromRows([][]float64{{1, 0}, {0, 1}})
jump := FromRows([][]float64{{1, 0}, {1, 0}})
reward := FromRows([][]float64{{1, 1}, {0, -0.1}})
m, _ := NewMDP([]*Matrix{stay, jump}, reward, 0.5)
res := m.ValueIteration(1e-12, 1000)
fmt.Printf("V = [%.2f %.2f], best action in state 1 = %d\n",
	res.Value[0], res.Value[1], res.Policy[1])

Output:
V = [2.00 0.90], best action in state 1 = 1

type Matrix 

type Matrix struct {
	// contains filtered or unexported fields
}

Matrix is a dense real matrix stored in row-major order. It is the basic linear-algebra container used throughout the optimalcontrol package. The zero value is not usable; construct matrices with NewMatrix, Zeros and friends.

func BlockMatrix 

func BlockMatrix(a, b, c, d *Matrix) *Matrix

BlockMatrix assembles the 2×2 block matrix [[a b];[c d]].

func CAREResidual 

func CAREResidual(a, b, q, r, x *Matrix) (*Matrix, error)

CAREResidual returns Aᵀ X + X A − X B R⁻¹ Bᵀ X + Q, the residual of the continuous algebraic Riccati equation for a candidate solution X.

func Cholesky 

func Cholesky(a *Matrix) (*Matrix, error)

Cholesky computes the lower-triangular Cholesky factor L with A = L Lᵀ for a symmetric positive-definite matrix. It returns an error if A is not SPD.

func ColumnVector 

func ColumnVector(v []float64) *Matrix

ColumnVector builds an n×1 matrix from v.

func ContinuousGain 

func ContinuousGain(b, r, p *Matrix) (*Matrix, error)

ContinuousGain returns the continuous LQR feedback gain K = R⁻¹ Bᵀ P.

func ControllabilityGramianContinuous 

func ControllabilityGramianContinuous(a, b *Matrix) (*Matrix, error)

ControllabilityGramianContinuous returns the infinite-horizon controllability Gramian Wc solving A Wc + Wc Aᵀ + B Bᵀ = 0 for a stable A.

func ControllabilityGramianDiscrete 

func ControllabilityGramianDiscrete(a, b *Matrix) (*Matrix, error)

ControllabilityGramianDiscrete returns the discrete controllability Gramian solving A Wc Aᵀ − Wc + B Bᵀ = 0 for a Schur-stable A.

func ControllabilityMatrix 

func ControllabilityMatrix(a, b *Matrix) *Matrix

ControllabilityMatrix returns the controllability matrix [B, AB, A²B, …, Aⁿ⁻¹B] for the pair (A, B).

func DAREResidual 

func DAREResidual(a, b, q, r, x *Matrix) (*Matrix, error)

DAREResidual returns the residual X − (Aᵀ X A − Aᵀ X B (R+BᵀXB)⁻¹ BᵀXA + Q) of the discrete algebraic Riccati equation for a candidate solution X.

func Diag 

func Diag(v []float64) *Matrix

Diag returns a square diagonal matrix whose diagonal is v.

func DiscreteGain 

func DiscreteGain(a, b, r, p *Matrix) (*Matrix, error)

DiscreteGain returns the discrete LQR feedback gain K = (R + Bᵀ P B)⁻¹ Bᵀ P A for a given Riccati solution P.

func DiscretizeZOH 

func DiscretizeZOH(a, b *Matrix, dt float64) (ad, bd *Matrix)

DiscretizeZOH returns the exact zero-order-hold discretization (Ad, Bd) of the continuous pair (A, B) for sample time dt, computed from the matrix exponential of the augmented block matrix [[A, B],[0, 0]].

func Eye 

func Eye(n int) *Matrix

Eye returns the n×n identity matrix.

func FromRows 

func FromRows(rows [][]float64) *Matrix

FromRows builds a matrix from a slice of equal-length rows.

func HStack 

func HStack(a, b *Matrix) *Matrix

HStack returns the horizontal concatenation [a b].

func HamiltonianMatrix 

func HamiltonianMatrix(a, b, q, r *Matrix) (*Matrix, error)

HamiltonianMatrix returns the 2n×2n Hamiltonian matrix 

H = [[ A,  -S ],
     [ -Q, -Aᵀ ]],  with  S = B R⁻¹ Bᵀ,

whose stable invariant subspace yields the stabilizing solution of the continuous-time algebraic Riccati equation.

func Inverse 

func Inverse(a *Matrix) (*Matrix, error)

Inverse returns the inverse of a square matrix.

func Kron 

func Kron(a, b *Matrix) *Matrix

Kron returns the Kronecker product a⊗b.

func LeastSquaresMatrix 

func LeastSquaresMatrix(a, b *Matrix) (*Matrix, error)

LeastSquaresMatrix solves A X = B in the least-squares sense column by column.

func LyapunovContinuousResidual 

func LyapunovContinuousResidual(a, q, x *Matrix) *Matrix

LyapunovContinuousResidual returns Aᵀ X + X A + Q, the residual of the continuous Lyapunov equation, useful for verifying a solution.

func LyapunovDiscreteResidual 

func LyapunovDiscreteResidual(a, q, x *Matrix) *Matrix

LyapunovDiscreteResidual returns Aᵀ X A − X + Q, the residual of the discrete Lyapunov equation.

func MatrixExp 

func MatrixExp(a *Matrix) *Matrix

MatrixExp returns the matrix exponential exp(A) using a scaling-and-squaring scheme with a truncated Taylor series. It is accurate for the moderate-sized matrices used in control applications.

func MatrixPow 

func MatrixPow(a *Matrix, p int) *Matrix

MatrixPow returns A raised to the non-negative integer power p.

func MatrixSign 

func MatrixSign(a *Matrix, maxIter int, tol float64) (*Matrix, error)

MatrixSign computes the matrix sign function of a square matrix with no purely imaginary eigenvalues, using the Newton iteration with determinantal scaling 

Z_{k+1} = ½ ( c_k Z_k + c_k⁻¹ Z_k⁻¹ ),   c_k = |det Z_k|^{-1/N}.

The returned matrix S satisfies S² = I and shares the eigenvectors of A, with eigenvalues ±1 according to the sign of the real part of A's eigenvalues.

func NewMatrix 

func NewMatrix(r, c int, data []float64) *Matrix

NewMatrix builds an r×c matrix from the supplied row-major data. The data slice is copied. It panics if len(data) != r*c.

func ObservabilityGramianContinuous 

func ObservabilityGramianContinuous(a, c *Matrix) (*Matrix, error)

ObservabilityGramianContinuous returns the infinite-horizon observability Gramian Wo solving Aᵀ Wo + Wo A + Cᵀ C = 0 for a stable A.

func ObservabilityGramianDiscrete 

func ObservabilityGramianDiscrete(a, c *Matrix) (*Matrix, error)

ObservabilityGramianDiscrete returns the discrete observability Gramian solving Aᵀ Wo A − Wo + Cᵀ C = 0 for a Schur-stable A.

func ObservabilityMatrix 

func ObservabilityMatrix(a, c *Matrix) *Matrix

ObservabilityMatrix returns the observability matrix [C; CA; CA²; …; CAⁿ⁻¹] for the pair (A, C).

func Ones 

func Ones(r, c int) *Matrix

Ones returns an r×c matrix whose entries are all one.

func RowVector 

func RowVector(v []float64) *Matrix

RowVector builds a 1×n matrix from v.

func SolveCARE 

func SolveCARE(a, b, q, r *Matrix) (*Matrix, error)

SolveCARE solves the continuous-time algebraic Riccati equation for its stabilizing solution. It is an alias for the robust matrix-sign method SolveCARESign.

func SolveCAREKleinman 

func SolveCAREKleinman(a, b, q, r, k0 *Matrix, maxIter int, tol float64) (*Matrix, error)

SolveCAREKleinman solves the continuous-time algebraic Riccati equation by Kleinman's Newton iteration starting from a stabilizing gain k0 (so that A − B k0 is Hurwitz). Each step solves a Lyapunov equation, giving quadratic convergence to the stabilizing solution.

func SolveCARESign 

func SolveCARESign(a, b, q, r *Matrix) (*Matrix, error)

SolveCARESign solves the continuous-time algebraic Riccati equation 

Aᵀ X + X A − X B R⁻¹ Bᵀ X + Q = 0

for the symmetric stabilizing solution X using the matrix-sign-function method applied to the Hamiltonian matrix. Q must be symmetric positive semidefinite and R symmetric positive definite; the pair (A, B) must be stabilizable and (A, Q) detectable.

func SolveDARE 

func SolveDARE(a, b, q, r *Matrix) (*Matrix, error)

SolveDARE solves the discrete-time algebraic Riccati equation for its stabilizing solution. It uses the iterative recursion SolveDAREIter with a generous iteration budget.

Example 
a := FromRows([][]float64{{1}})
b := FromRows([][]float64{{1}})
q := FromRows([][]float64{{1}})
r := FromRows([][]float64{{1}})
p, _ := SolveDARE(a, b, q, r)
fmt.Printf("P = %.4f\n", p.At(0, 0)) // golden ratio

Output:
P = 1.6180

func SolveDAREIter 

func SolveDAREIter(a, b, q, r *Matrix, maxIter int, tol float64) (*Matrix, error)

SolveDAREIter solves the discrete-time algebraic Riccati equation 

X = Aᵀ X A − Aᵀ X B (R + Bᵀ X B)⁻¹ Bᵀ X A + Q

by fixed-point iteration of the Riccati recursion, starting from X₀ = Q.

func SolveLyapunovContinuous 

func SolveLyapunovContinuous(a, q *Matrix) (*Matrix, error)

SolveLyapunovContinuous solves the continuous-time Lyapunov equation Aᵀ X + X A + Q = 0 for the symmetric matrix X.

func SolveLyapunovDiscrete 

func SolveLyapunovDiscrete(a, q *Matrix) (*Matrix, error)

SolveLyapunovDiscrete solves the discrete-time (Stein) Lyapunov equation Aᵀ X A − X + Q = 0 for the symmetric matrix X.

func SolveMatrix 

func SolveMatrix(a, b *Matrix) (*Matrix, error)

SolveMatrix solves A X = B for the matrix X.

func SolveSylvester 

func SolveSylvester(a, b, c *Matrix) (*Matrix, error)

SolveSylvester solves the Sylvester equation A X + X B = C for X, where A is m×m, B is n×n and C is m×n. It uses the Kronecker-product formulation (I⊗A + Bᵀ⊗I) vec(X) = vec(C) and a dense linear solve, which is robust for the small systems arising in control design.

func SymEigen 

func SymEigen(a *Matrix) (w []float64, v *Matrix, err error)

SymEigen computes eigenvalues and orthonormal eigenvectors of a symmetric matrix using the cyclic Jacobi method. The eigenvalues are returned as a slice and the eigenvectors as the columns of V, so that A V = V diag(w).

func Unvec 

func Unvec(v []float64, r, c int) *Matrix

Unvec reshapes a column-major vector into an r×c matrix.

func VStack 

func VStack(a, b *Matrix) *Matrix

VStack returns the vertical concatenation [a; b].

func Zeros 

func Zeros(r, c int) *Matrix

Zeros returns an r×c matrix of zeros.

func (*Matrix) Add 

func (m *Matrix) Add(i, j int, v float64)

Add accumulates v into the element at row i, column j.

func (*Matrix) ApproxEqual 

func (m *Matrix) ApproxEqual(b *Matrix, tol float64) bool

ApproxEqual reports whether a and b have the same shape and entries within absolute tolerance tol.

func (*Matrix) At 

func (m *Matrix) At(i, j int) float64

At returns the element at row i, column j.

func (*Matrix) Clone 

func (m *Matrix) Clone() *Matrix

Clone returns a deep copy of the matrix.

func (*Matrix) Col 

func (m *Matrix) Col(j int) []float64

Col returns a copy of column j as a slice.

func (*Matrix) Cols 

func (m *Matrix) Cols() int

Cols returns the number of columns.

func (*Matrix) Data 

func (m *Matrix) Data() []float64

Data returns a copy of the underlying row-major data.

func (*Matrix) Equal 

func (m *Matrix) Equal(b *Matrix) bool

Equal reports whether a and b have the same shape and identical entries.

func (*Matrix) FrobeniusNorm 

func (m *Matrix) FrobeniusNorm() float64

FrobeniusNorm returns the Frobenius norm of the matrix.

func (*Matrix) InfNorm 

func (m *Matrix) InfNorm() float64

InfNorm returns the maximum absolute row sum of the matrix.

func (*Matrix) IsSquare 

func (m *Matrix) IsSquare() bool

IsSquare reports whether the matrix is square.

func (*Matrix) IsSymmetric 

func (m *Matrix) IsSymmetric(tol float64) bool

IsSymmetric reports whether the matrix is square and symmetric within tol.

func (*Matrix) MaxAbs 

func (m *Matrix) MaxAbs() float64

MaxAbs returns the largest absolute entry of the matrix.

func (*Matrix) Minus 

func (m *Matrix) Minus(b *Matrix) *Matrix

Minus returns the difference a-b.

func (*Matrix) Mul 

func (m *Matrix) Mul(b *Matrix) *Matrix

Mul returns the matrix product a·b.

func (*Matrix) MulVec 

func (m *Matrix) MulVec(v []float64) []float64

MulVec returns the matrix-vector product a·v.

func (*Matrix) Neg 

func (m *Matrix) Neg() *Matrix

Neg returns the additive inverse of the matrix.

func (*Matrix) OneNorm 

func (m *Matrix) OneNorm() float64

OneNorm returns the maximum absolute column sum of the matrix.

func (*Matrix) Plus 

func (m *Matrix) Plus(b *Matrix) *Matrix

Plus returns the sum a+b.

func (*Matrix) Row 

func (m *Matrix) Row(i int) []float64

Row returns a copy of row i as a slice.

func (*Matrix) Rows 

func (m *Matrix) Rows() int

Rows returns the number of rows.

func (*Matrix) Scale 

func (m *Matrix) Scale(s float64) *Matrix

Scale returns the matrix scaled by s.

func (*Matrix) Set 

func (m *Matrix) Set(i, j int, v float64)

Set assigns v to the element at row i, column j.

func (*Matrix) SetBlock 

func (m *Matrix) SetBlock(r0, c0 int, b *Matrix)

SetBlock copies b into the block of m with top-left corner at (r0, c0).

func (*Matrix) SetCol 

func (m *Matrix) SetCol(j int, v []float64)

SetCol overwrites column j with v.

func (*Matrix) SetRow 

func (m *Matrix) SetRow(i int, v []float64)

SetRow overwrites row i with v.

func (*Matrix) String 

func (m *Matrix) String() string

String renders the matrix for debugging.

func (*Matrix) Submatrix 

func (m *Matrix) Submatrix(r0, r1, c0, c1 int) *Matrix

Submatrix returns the r0..r1-1 × c0..c1-1 block of the matrix.

func (*Matrix) Symmetrize 

func (m *Matrix) Symmetrize() *Matrix

Symmetrize returns (A + Aᵀ)/2, the symmetric part of a square matrix.

func (*Matrix) Trace 

func (m *Matrix) Trace() float64

Trace returns the sum of the diagonal entries of a square matrix.

func (*Matrix) Transpose 

func (m *Matrix) Transpose() *Matrix

Transpose returns the transpose of the matrix.

func (*Matrix) Vec 

func (m *Matrix) Vec() []float64

Vec returns the column-major vectorization of the matrix (columns stacked).

type PolicyIterationResult 

type PolicyIterationResult struct {
	// Policy is the optimal deterministic policy.
	Policy []int
	// Value is the value function of the optimal policy.
	Value []float64
	// Iterations is the number of policy-improvement steps performed.
	Iterations int
}

PolicyIterationResult holds the outcome of policy iteration.

type TPBVPSolution 

type TPBVPSolution struct {
	// Times holds the ascending sample times from 0 to T.
	Times []float64
	// X holds the state at each sample time.
	X [][]float64
	// P holds the costate (adjoint) at each sample time.
	P [][]float64
	// U holds the optimal control at each sample time.
	U [][]float64
}

TPBVPSolution holds the sampled solution of a two-point boundary-value problem arising from Pontryagin's minimum principle: the state, costate and optimal control trajectories on a time grid.

func IndirectShooting 

func IndirectShooting(
	x0 []float64,
	f func(x, u []float64) []float64,
	g func(x, p, u []float64) []float64,
	uOpt func(x, p []float64) []float64,
	termGrad func(x []float64) []float64,
	tFinal float64, steps, maxIter int, tol float64,
) ([]float64, *TPBVPSolution, error)

IndirectShooting solves a fixed-final-time two-point boundary-value problem by single shooting on the initial costate. The problem is specified by the augmented dynamics 

x' = f(x, u),   p' = g(x, p, u),   u = uOpt(x, p),

with x(0) = x0 and terminal condition p(T) = termGrad(x(T)). Newton's method with a finite-difference Jacobian is applied to the shooting residual. It returns the initial costate and the sampled trajectory.

func PontryaginLQ 

func PontryaginLQ(a, b, q, r, qf *Matrix, x0 []float64, tFinal float64, steps int) (*TPBVPSolution, error)

PontryaginLQ solves the fixed-final-time linear-quadratic optimal-control problem 

min  ½ x(T)ᵀ Qf x(T) + ½ ∫₀ᵀ (xᵀ Q x + uᵀ R u) dt
s.t. x' = A x + B u,  x(0) = x0,

via Pontryagin's minimum principle. It forms the Hamiltonian two-point boundary-value problem, propagates it with the Hamiltonian state-transition matrix, solves the terminal transversality condition p(T) = Qf x(T) for the initial costate, and returns the sampled state, costate and control.

type ValueIterationResult 

type ValueIterationResult struct {
	// Value is the (approximately) optimal value function.
	Value []float64
	// Policy is the greedy policy with respect to Value.
	Policy []int
	// Iterations is the number of sweeps performed.
	Iterations int
	// Converged reports whether the max-norm residual fell below the tolerance.
	Converged bool
}

ValueIterationResult holds the outcome of value iteration.

Source Files

View all Source files

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL
go.dev uses cookies from Google to deliver and enhance the quality of its services and to analyze traffic. Learn more.