SOTN-001

OFC Control Loop Mathematical Description#

Abstract

This technical note provides a comprehensive mathematical description of the Optical Feedback Control (OFC) system used in the Vera C. Rubin Observatory’s Active Optics System (AOS), as implemented in ts_ofc and ts_mtaos.

The document begins with the forward optical model, explaining how the sensitivity matrix relates telescope degrees of freedom (DOFs) to wavefront Zernike coefficients, and how DOF normalization ensures balanced contributions during state estimation.

The core of the note addresses the state estimation inverse problem: recovering DOF states from noisy wavefront measurements. Two complementary regularization approaches are presented — SVD truncation and noise-weighted least squares — along with a concrete numerical example showing matrix dimensions for the typical 76-measurement, 22-DOF configuration.

A dedicated section defines v-modes (optical eigenmodes from the SVD), explains their physical interpretation as coordinated DOF adjustments that produce characteristic aberration patterns, and describes coordinate transformations between DOF space and v-mode space.

Finally, the note covers PID control in both DOF space and v-mode space, showing how v-mode gains enable tuning by optical importance, and how corrections are aggregated and dispatched to the telescope subsystems.

This document serves as a reference for developers working on AOS software and operators seeking to understand the control loop mathematics.

V-Modes and the OFC Control Loop#

This chapter provides a mathematical description of the Optical Feedback Control (OFC) pipeline implemented in ts_ofc and ts_mtaos. It covers the forward optical model, the state estimation inverse problem, the definition and physical interpretation of v-modes, and how PID gains are applied in both DOF space and v-mode space.

The forward optical model#

The relationship between the telescope’s degrees of freedom (DOFs) and the resulting wavefront error across the focal plane is captured by the sensitivity matrix. For each sensor position on the focal plane, the sensitivity matrix encodes how a unit change in each DOF affects each Zernike coefficient of the wavefront.

If we stack the Zernike coefficients from all sensors into a single vector, the forward model is:

\[\mathbf{y} = \mathbf{A}\,\mathbf{x}\]

where:

  • \(\mathbf{y} \in \mathbb{R}^{m}\) is the wavefront error vector (Zernike coefficients across all sensors), with \(m = n_\text{sensors} \times n_\text{Zernikes}\). For four corner wavefront sensors and Zernikes Z4–Z22 (19 coefficients), \(m = 4 \times 19 = 76\).

  • \(\mathbf{A} \in \mathbb{R}^{m \times n}\) is the sensitivity matrix, evaluated at the sensor field positions and the current camera rotation angle.

  • \(\mathbf{x} \in \mathbb{R}^{n}\) is the DOF state vector, with \(n = n_\text{used DOFs}\). A typical configuration uses \(n = 22\) DOFs: 5 M2 hexapod, 5 camera hexapod, 7 M1M3 bending modes, and 5 M2 bending modes.

The sensitivity matrix \(\mathbf{A}\) is constructed from a double Zernike expansion of the optical model and is evaluated at each sensor’s field angle. It depends on the camera rotation angle because the Zernike basis rotates with the instrument.

DOF normalization#

The physical DOFs span different scales: hexapod translations are in micrometers, hexapod rotations in degrees, and bending mode coefficients are dimensionless with varying magnitudes. To give each DOF comparable influence during the inversion, we introduce a diagonal normalization matrix:

\[\mathbf{N} = \operatorname{diag}\!\bigl(w_1, w_2, \ldots, w_n\bigr)\]

where the weights \(w_i\) are read from the normalization_weights configuration file (for details see Megias Homar, G. et al. 2024, Eqs. 9–11).

The normalized sensitivity matrix is:

\[\tilde{\mathbf{A}} = \mathbf{A}\,\mathbf{N}\]

so that the forward model in normalized DOF coordinates \(\tilde{\mathbf{x}} = \mathbf{N}^{-1}\mathbf{x}\) becomes:

\[\mathbf{y} = \tilde{\mathbf{A}}\,\tilde{\mathbf{x}}\]

State estimation: the inverse problem#

Given a measured wavefront error \(\mathbf{y}_\text{meas}\), we want to estimate the DOF state \(\hat{\mathbf{x}}\). This is an overdetermined linear inverse problem (\(m > n\)), which can be solved via regularized pseudo-inversion. There are two main approaches: SVD truncation and noise-weighted least squares. The OFC combines both for optimal state estimation.

Problem setup#

The residual wavefront error used as input is:

\[\mathbf{y} = \mathbf{y}_\text{meas} - \mathbf{y}_\text{intrinsic} - \mathbf{y}_2\]

These subtractions operate in measurement space (Zernike coefficients at each sensor, \(\mathbf{y} \in \mathbb{R}^m\)):

  • \(\mathbf{y}_\text{intrinsic}\) contains the design (intrinsic) Zernike coefficients for the current filter and sensor positions. These represent the wavefront error that the optical design produces even when the optics are in perfect alignment (\(\mathbf{x} = \mathbf{0}\)). When subtract_intrinsics = True (the default for the main closed-loop path), these are subtracted so the residual \(\mathbf{y}\) reflects only the misalignment component.

  • \(\mathbf{y}_2\) is a static per-sensor correction in measurement space. In the deployed configuration (lsst_y2.yaml), only the Z4 (focus) entries for the corner wavefront sensors are non-zero (e.g., R00: -0.139 um, R44: -0.145 um). These values were determined from on-sky data as the focus offsets at which the best image quality is achieved.

The purpose of these subtractions is to ensure that when the telescope is at the desired alignment state (\(\mathbf{x} = \mathbf{0}\) in DOF space), the residual \(\mathbf{y} \approx \mathbf{0}\). The state estimator then inverts \(\mathbf{y}\) to produce the DOF estimate \(\hat{\mathbf{x}}\). The PID controller operates in DOF space (\(\mathbf{x} \in \mathbb{R}^n\)) with a setpoint of zero (see PID control below). So the full chain is:

  1. Subtract known contributions in y-space: \(\mathbf{y} = \mathbf{y}_\text{meas} - \mathbf{y}_\text{intrinsic} - \mathbf{y}_2\)

  2. Estimate the DOF state: \(\hat{\mathbf{x}} = f(\mathbf{y})\)

  3. PID drives \(\hat{\mathbf{x}} \to \mathbf{0}\) (the setpoint in x-space)

The y-space subtractions and the x-space setpoint work together: the subtractions define what “zero misalignment” looks like in measurement space, and the PID targets zero in DOF space. If the y2 values were wrong (e.g., the focus offsets didn’t correspond to best image quality), the system would converge to a suboptimal state even though the PID error would be zero.

We seek to solve:

\[\mathbf{y} = \tilde{\mathbf{A}}\,\tilde{\mathbf{x}} + \mathbf{n}\]

where \(\tilde{\mathbf{A}}\) is the normalized sensitivity matrix, \(\tilde{\mathbf{x}}\) is the normalized DOF state, and \(\mathbf{n}\) is measurement noise with covariance \(\mathbf{C}_n = \mathbb{E}[\mathbf{n}\mathbf{n}^T]\).

Here \(\mathbb{E}[\cdot]\) denotes the expectation operator (statistical mean), and \(\mathbf{n}\mathbf{n}^T\) is the outer product of the noise vector with itself. The covariance matrix \(\mathbf{C}_n\) therefore contains:

  • Diagonal entries \((\mathbf{C}_n)_{ii} = \mathbb{E}[n_i^2] = \sigma_i^2\): the variance of the \(i\)-th Zernike measurement.

  • Off-diagonal entries \((\mathbf{C}_n)_{ij} = \mathbb{E}[n_i n_j]\): the covariance between Zernikes \(i\) and \(j\).

In practice, \(\mathbf{C}_n\) is estimated from stability data by computing the sample covariance of the residual Zernikes across \(N\) exposures:

\[\hat{\mathbf{C}}_n = \frac{1}{N-1} \sum_{k=1}^{N} (\mathbf{y}_k - \bar{\mathbf{y}}) (\mathbf{y}_k - \bar{\mathbf{y}})^T\]

where \(\mathbf{y}_k\) are the measured Zernike vectors and \(\bar{\mathbf{y}}\) is their mean. This matrix tells the state estimator which Zernikes and corners are noisy (large diagonal entries) and which are correlated (large off-diagonal entries), enabling appropriate down-weighting during state estimation. See SITCOMTN-177 for details on how this covariance is estimated from on-sky data.

Approach 1: SVD with truncation#

The standard least-squares solution uses the pseudo-inverse. Start by computing the singular value decomposition (SVD):

\[\tilde{\mathbf{A}} = \mathbf{U}\,\operatorname{diag}(\boldsymbol{\Sigma})\,\mathbf{V}_h\]

In the full SVD of an \(m \times n\) matrix with \(m > n\):

  • \(\mathbf{U} \in \mathbb{R}^{m \times m}\) is an orthogonal matrix (left singular vectors),

  • \(\boldsymbol{\Sigma} \in \mathbb{R}^{n}\) is a vector of \(n\) singular values \((\sigma_1, \ldots, \sigma_n)\), and \(\operatorname{diag}(\boldsymbol{\Sigma})\) is the \(m \times n\) rectangular matrix with these values on the main diagonal,

  • \(\mathbf{V}_h \in \mathbb{R}^{n \times n}\) is an orthogonal matrix (right singular vectors).

The singular values are ordered \(\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_n \geq 0\).

In practice, the thin (economy) SVD is used, which discards the extra columns of \(\mathbf{U}\) corresponding to zero rows of the diagonal matrix:

  • \(\mathbf{U} \in \mathbb{R}^{m \times n}\) — only the first \(n\) columns,

  • \(\boldsymbol{\Sigma} \in \mathbb{R}^{n}\) — vector of \(n\) singular values \((\sigma_1, \sigma_2, \ldots, \sigma_n)\),

  • \(\mathbf{V}_h \in \mathbb{R}^{n \times n}\) — unchanged.

In this case, \(\operatorname{diag}(\boldsymbol{\Sigma})\) is an \(n \times n\) square diagonal matrix.

The pseudo-inverse is:

\[\tilde{\mathbf{A}}^{+} = \mathbf{V}_h^T\,\operatorname{diag}(1/\boldsymbol{\Sigma})\,\mathbf{U}^T\]

where \(\operatorname{diag}(1/\boldsymbol{\Sigma}) = \operatorname{diag}(1/\sigma_1, \ldots, 1/\sigma_n)\).

Truncation regularization. Small singular values amplify noise because \(1/\sigma_i \to \infty\) as \(\sigma_i \to 0\). To regularize, truncate at index \(k\) by keeping only the first \(k\) singular values:

\[\tilde{\mathbf{A}}_k^{+} = \mathbf{V}_{h,k}^T\,\operatorname{diag}(1/\boldsymbol{\Sigma}_k)\,\mathbf{U}_k^T\]

where \(\mathbf{U}_k\) contains the first \(k\) columns of \(\mathbf{U}\), \(\mathbf{V}_{h,k}\) contains the first \(k\) rows of \(\mathbf{V}_h\), and \(\boldsymbol{\Sigma}_k\) contains the first \(k\) elements of \(\boldsymbol{\Sigma}\).

The truncated least-squares solution is:

\[\tilde{\mathbf{x}}_\text{trunc} = \tilde{\mathbf{A}}_k^{+}\,\mathbf{y}\]

This discards the \(n - k\) least-constrained DOF combinations, retaining only directions in DOF space that produce measurable wavefront changes above the noise floor.

Concrete example. Consider the typical OFC configuration with \(m = 76\) measurements (4 corners × 19 Zernikes) and \(n = 22\) DOFs. The thin SVD of \(\tilde{\mathbf{A}} \in \mathbb{R}^{76 \times 22}\) produces:

  • \(\mathbf{U} \in \mathbb{R}^{76 \times 22}\) — 22 left singular vectors (columns)

  • \(\boldsymbol{\Sigma} \in \mathbb{R}^{22}\) — vector of 22 singular values \((\sigma_1, \ldots, \sigma_{22})\)

  • \(\mathbf{V}_h \in \mathbb{R}^{22 \times 22}\) — 22 right singular vectors (the v-modes)

With truncation_index = 12, we keep only the first \(k = 12\) modes:

  • \(\mathbf{U}_{12} \in \mathbb{R}^{76 \times 12}\) — first 12 columns of U

  • \(\boldsymbol{\Sigma}_{12} \in \mathbb{R}^{12}\) — first 12 singular values \((\sigma_1, \ldots, \sigma_{12})\)

  • \(\mathbf{V}_{h,12} \in \mathbb{R}^{12 \times 22}\) — first 12 rows of Vₕ

The truncated pseudo-inverse is then:

\[\tilde{\mathbf{A}}_{12}^{+} = \underbrace{\mathbf{V}_{h,12}^T}_{22 \times 12} \underbrace{\operatorname{diag}(1/\boldsymbol{\Sigma}_{12})}_{12 \times 12} \underbrace{\mathbf{U}_{12}^T}_{12 \times 76} \;\in\; \mathbb{R}^{22 \times 76}\]

The final state estimate is:

\[\underbrace{\tilde{\mathbf{x}}}_{22 \times 1} = \underbrace{\tilde{\mathbf{A}}_{12}^{+}}_{22 \times 76} \cdot \underbrace{\mathbf{y}}_{76 \times 1}\]

The output \(\tilde{\mathbf{x}}\) has 22 components (one per DOF), but only 12 independent directions are actually constrained by the data. Internally, the matrix multiplication proceeds right-to-left:

  1. \(\mathbf{U}_{12}^T \cdot \mathbf{y}\) projects the 76 measurements onto 12 left singular vectors, producing a 12-element vector.

  2. \(\operatorname{diag}(1/\boldsymbol{\Sigma}_{12})\) scales each component by the inverse of its singular value.

  3. \(\mathbf{V}_{h,12}^T\) expands this 12-element vector back to 22 DOFs using the first 12 v-modes.

The result lies entirely in the subspace spanned by v-modes 0–11. The coefficients of v-modes 12–21 are implicitly zero, preventing noise amplification from small singular values \(\sigma_{13}\) through \(\sigma_{22}\).

Approach 2: Noise-weighted pseudo-inverse#

When the measurement noise covariance \(\mathbf{C}_n\) is known and non-uniform, we want to minimize the weighted least-squares objective:

\[\min_{\tilde{\mathbf{x}}} \; \bigl(\mathbf{y} - \tilde{\mathbf{A}}\,\tilde{\mathbf{x}}\bigr)^T \mathbf{C}_n^{-1} \bigl(\mathbf{y} - \tilde{\mathbf{A}}\,\tilde{\mathbf{x}}\bigr)\]

The solution is:

\[\tilde{\mathbf{x}}_\text{wls} = \bigl(\tilde{\mathbf{A}}^T \mathbf{C}_n^{-1} \tilde{\mathbf{A}}\bigr)^{-1} \tilde{\mathbf{A}}^T \mathbf{C}_n^{-1} \mathbf{y}\]

Equivalently, whiten the problem by defining:

\[\bar{\mathbf{y}} = \mathbf{C}_n^{-1/2}\,\mathbf{y}, \qquad \bar{\mathbf{A}} = \mathbf{C}_n^{-1/2}\,\tilde{\mathbf{A}}\]

Then solve the standard (unweighted) least-squares:

\[\tilde{\mathbf{x}}_\text{wls} = \bar{\mathbf{A}}^{+}\,\bar{\mathbf{y}}\]

The noise covariance down-weights Zernike modes and sensors with high measurement uncertainty. See SITCOMTN-177 for details on how \(\mathbf{C}_n\) is estimated from on-sky stability data.

Combined approach (as implemented in ts_ofc)#

The StateEstimator in ts_ofc combines both approaches:

Step 1 — Whiten by noise covariance:

\[\bar{\mathbf{y}} = \mathbf{C}_n^{-1/2}\,\mathbf{y}, \qquad \bar{\mathbf{A}} = \mathbf{C}_n^{-1/2}\,\tilde{\mathbf{A}}\]

Step 2 — Truncated pseudo-inverse of whitened system:

\[\tilde{\mathbf{x}} = \bar{\mathbf{A}}_k^{+}\,\bar{\mathbf{y}}\]

where \(\bar{\mathbf{A}}_k^{+}\) is the truncated pseudo-inverse of the whitened sensitivity matrix (SVD with \(k\) retained singular values).

Step 3 — De-normalize to physical DOFs:

\[\hat{\mathbf{x}} = \mathbf{N}\,\tilde{\mathbf{x}}\]

In compact form:

\[\hat{\mathbf{x}} = \mathbf{N} \cdot \bigl(\mathbf{C}_n^{-1/2}\,\tilde{\mathbf{A}}\bigr)_k^{+} \cdot \mathbf{C}_n^{-1/2}\,\mathbf{y}\]

For implementation details, see Equation (10) in arXiv:2406.04656 and Martínez-Galarza et al. 2024.

Definition of v-modes#

The v-modes are defined as the right singular vectors of the normalized sensitivity matrix. From the thin SVD:

\[\tilde{\mathbf{A}} = \mathbf{U}\,\operatorname{diag}(\boldsymbol{\Sigma})\,\mathbf{V}_h\]

the rows of \(\mathbf{V}_h\) form an orthonormal basis for the normalized DOF space. Row \(i\) of \(\mathbf{V}_h\) is v-mode \(i\).

Each v-mode is a linear combination of all (normalized) DOFs. Concretely, v-mode \(i\) is the vector \(\mathbf{v}_i = (\mathbf{V}_h)_{i,:} \in \mathbb{R}^{n}\), which specifies a coordinated adjustment of hexapod positions and bending modes.

The associated singular value \(\sigma_i\) quantifies the optical leverage of v-mode \(i\): the amount of wavefront change (in Zernike amplitude) produced per unit of that DOF combination. Because the singular values are sorted in decreasing order, v-mode 0 is the most optically important mode and v-mode \(n{-}1\) is the least.

The number of v-modes equals the number of used DOFs. With \(n = 22\) used DOFs, the SVD of the \((76 \times 22)\) matrix produces \(\mathbf{V}_h \in \mathbb{R}^{22 \times 22}\), yielding exactly 22 v-modes.

Coordinate transformations#

DOFs to v-modes. Given a DOF vector \(\mathbf{x}\), the v-mode representation is:

\[\mathbf{v} = \mathbf{V}_h\,\mathbf{N}^{-1}\,\mathbf{x}\]

V-modes to DOFs. Given a v-mode vector \(\mathbf{v}\), the DOF representation is:

\[\mathbf{x} = \mathbf{N}\,\mathbf{V}_h^{T}\,\mathbf{v}\]

These transformations are exact and lossless because \(\mathbf{V}_h^{T}\mathbf{V}_h = \mathbf{I}\).

In the code, these correspond to StateEstimator.get_vmodes_from_dofs() and StateEstimator.get_dofs_from_vmodes() in state_estimator.py.

Physical interpretation of v-modes#

Although each v-mode is a linear combination of all DOFs, each v-mode produces a characteristic wavefront pattern described by the corresponding left singular vector \(\mathbf{u}_i = \mathbf{U}_{:,i}\). This vector contains the Zernike coefficients (across all sensors) that result from commanding v-mode \(i\).

If \(\mathbf{u}_i\) has most of its energy in the Zernike coefficients corresponding to a particular aberration — say, focus (Z4) — then v-mode \(i\) is effectively “the focus mode”: the coordinated DOF adjustment that most efficiently changes the telescope’s focus. Similarly, if \(\mathbf{u}_j\) is dominated by spherical aberration (Z11), then v-mode \(j\) is “the spherical mode.”

This is why operators associate specific v-modes with specific aberrations:

  • V-mode \(i\) is named after the dominant Zernike in \(\mathbf{u}_i\).

  • The association is approximate — each v-mode typically produces a dominant Zernike plus smaller contributions to other Zernikes.

  • The SVD discovers these “optical eigenmodes” automatically from the sensitivity matrix; they are not prescribed by the user.

Truncation and v-mode importance#

The singular values \(\sigma_i\) determine which v-modes are well-constrained by the wavefront measurements:

  • V-modes \(0\) through \(k{-}1\) (large \(\sigma_i\)): These DOF combinations produce large, easily measured wavefront changes. The state estimator can reliably reconstruct their amplitudes.

  • V-modes \(k\) through \(n{-}1\) (small \(\sigma_i\)): These combinations produce tiny wavefront changes that are buried in measurement noise. The truncation at index \(k\) discards these during state estimation.

With truncation_index = 12 and \(n = 22\) DOFs, the first 12 v-modes are retained and the remaining 10 are discarded during state estimation. This does not reduce the dimension of the correction — the PID controller still operates on all 22 v-modes (or DOFs) — but the state estimate for v-modes 12–21 is effectively zero because their singular values were zeroed in the truncated sensitivity matrix.

PID control#

The PID controller computes a correction \(\mathbf{u}_k\) at each iteration \(k\) based on the estimated state. The correction is applied with negative feedback: the telescope moves in the opposite direction of the estimated error.

The overall flow in OFC.calculate_corrections() is:

  1. Estimate the optical state \(\hat{\mathbf{x}}_k\) via the state estimator.

  2. Pass \(\hat{\mathbf{x}}_k\) to the PID controller, which computes an internal control signal.

  3. Negate the control signal: uk = -controller.control_step(...)

  4. Aggregate the correction into the telescope state.

The negation in step 3 ensures that if the estimated state indicates a positive defocus, the correction drives it back toward zero.

The setpoint: what the PID targets#

The PID controller drives the estimated state toward a setpoint vector \(\mathbf{s} \in \mathbb{R}^{50}\), which specifies the desired value for each of the 50 possible degrees of freedom. In the default configuration (pid_controller.yaml), all setpoint entries are zero:

\[\mathbf{s} = \mathbf{0}\]

This means the controller targets a state where all DOF deviations are zero — i.e., the telescope optics are in their nominal alignment. A non-zero setpoint would allow the controller to target a deliberate offset (e.g., intentional defocus for engineering purposes).

The effective target of the closed loop is not simply “zero physical DOFs.” The PID setpoint lives in DOF space (\(\mathbf{x}\)), but what “zero” means depends on the measurement preprocessing (see Problem setup above). Before state estimation, the intrinsic Zernikes and y2 correction are subtracted in measurement space (\(\mathbf{y}\)). The y2 correction encodes per-sensor focus offsets determined from on-sky data.

The system converges when \(\hat{\mathbf{x}} = \mathbf{0}\), which corresponds to the measurement residual \(\mathbf{y} \approx \mathbf{0}\), which in turn means \(\mathbf{y}_\text{meas} \approx \mathbf{y}_\text{intrinsic} + \mathbf{y}_2\). So the telescope is driven to the state where its measured wavefront matches the intrinsic design prediction plus the empirical y2 offsets — i.e., the state that produces the best image quality as determined from on-sky calibration.

Only the entries corresponding to the active DOFs (dof_idx) are used.

PID formula (as implemented)#

At each iteration \(k\), the controller receives the estimated state \(\hat{\mathbf{x}}_k\) (either in DOF or v-mode space) and computes:

Error signal:

\[\mathbf{e}_k = \mathbf{s}[\text{dof\_idx}] - \hat{\mathbf{x}}_k\]

Integral with anti-windup:

\[\begin{split}\mathbf{I}_k &= \mathbf{I}_{k-1} + \mathbf{e}_k \\ \mathbf{I}_k &\leftarrow \operatorname{clip}\bigl( \mathbf{I}_k,\; -\mathbf{I}_\text{max}[\text{dof\_idx}],\; +\mathbf{I}_\text{max}[\text{dof\_idx}] \bigr)\end{split}\]

The max_integral vector (default: 1.0 for all DOFs) prevents the integral from growing without bound during sustained errors (anti-windup protection).

Derivative with low-pass filter:

\[\begin{split}\mathbf{d}_k^\text{raw} &= \mathbf{e}_k - \mathbf{e}_{k-1} \\ \mathbf{d}_k &= \alpha\,\mathbf{d}_k^\text{raw} + (1 - \alpha)\,\mathbf{d}_{k-1}\end{split}\]

where \(\alpha\) is the derivative_filter_coeff. When \(\alpha = 1\) (the default), no filtering is applied and the derivative is simply the difference between consecutive errors. When \(\alpha < 1\), the derivative is smoothed with an exponential moving average to reduce sensitivity to noise.

Control output:

\[\mathbf{u}_k^\text{pid} = \mathbf{K}_p\,\mathbf{e}_k + \mathbf{K}_i\,\mathbf{I}_k + \mathbf{K}_d\,\mathbf{d}_k\]

where:

  • \(\mathbf{K}_p = \operatorname{diag}\!\bigl(k_{p,1}, \ldots, k_{p,n}\bigr)\) is the proportional gain (default: 0.7 for all DOFs),

  • \(\mathbf{K}_i = \operatorname{diag}\!\bigl(k_{i,1}, \ldots, k_{i,n}\bigr)\) is the integral gain (default: 0.0 — integral action disabled),

  • \(\mathbf{K}_d = \operatorname{diag}\!\bigl(k_{d,1}, \ldots, k_{d,n}\bigr)\) is the derivative gain (default: 0.1).

The gains can be a scalar (applied uniformly to all DOFs) or an array (per-DOF tuning). They are indexed by dof_idx to select only the active DOFs.

The final correction applied to the telescope is:

\[\mathbf{u}_k = -\mathbf{u}_k^\text{pid}\]

This negation, applied in OFC.calculate_corrections(), converts the PID output (which points from state toward setpoint) into a physical correction (which drives the optics back toward the setpoint).

Default configuration#

The default PID parameters from pid_controller.yaml:

Parameter

Default

Description

kp

0.7

Proportional gain (scalar or per-DOF array)

ki

0.0

Integral gain (disabled by default)

kd

0.1

Derivative gain

derivative_filter_coeff

1.0

Low-pass filter coefficient (1.0 = no filtering)

setpoint

all zeros (50 elements)

Target state for each DOF

max_integral

1.0 (all DOFs)

Anti-windup clipping limit

Physical meaning of gains: \(k_p = 0.7\) means “correct 70% of the estimated error per iteration.” \(k_d = 0.1\) provides damping that reduces overshoot by opposing rapid changes in the error. With \(k_i = 0\), there is no integral action by default — steady-state errors are not automatically eliminated across iterations.

Controller state and history#

The controller maintains internal state between iterations:

  • aggregated_state (\(\mathbf{x}_\text{agg}\)): the cumulative correction applied to the telescope (initialized from dof_state0).

  • integral (\(\mathbf{I}\)): accumulated integral of the error.

  • previous_error (\(\mathbf{e}_{k-1}\)): error from the previous step.

  • filtered_derivative (\(\mathbf{d}_{k-1}\)): filtered derivative from the previous step.

The reset_history() method:

  1. Saves the current aggregated state as the new reference: \(\mathbf{x}_0 \leftarrow \mathbf{x}_\text{agg}\).

  2. Zeros the integral, previous error, and filtered derivative.

This is called when the closed loop is restarted with a new comp_dof_idx configuration, effectively starting the PID fresh from the current telescope state.

PID in DOF space#

When control_vmodes = False, the state estimator returns the optical state directly in DOF space:

\[\hat{\mathbf{x}}_k = \mathbf{N} \cdot \bigl(\mathbf{C}_n^{-1/2}\,\tilde{\mathbf{A}}\bigr)_k^{+} \cdot \mathbf{C}_n^{-1/2}\,\mathbf{y}\]

The PID operates on this 22-dimensional (or however many DOFs are active) state vector. Each gain entry corresponds directly to a physical DOF.

Physical meaning of DOF gains: each diagonal entry of \(\mathbf{K}_p\) controls the aggressiveness of correction for one physical DOF. For example, \(k_{p,0} = 0.7\) means “correct 70% of the estimated M2 hexapod dZ error per iteration.”

PID in v-mode space#

When control_vmodes = True, the state estimator returns the optical state directly in v-mode space. Internally, the state estimator solves the least-squares problem and returns the v-mode coefficients (zero-padded to the full DOF dimension):

\[\hat{\mathbf{v}}_k = \text{pad}_{n}\bigl(\mathbf{v}_\text{state}\bigr)\]

where \(\mathbf{v}_\text{state}\) has dimension equal to the number of retained modes (truncation index), padded with zeros for the discarded modes.

The PID then operates on this v-mode vector using the same calculate_pid_step() method — the controller does not know whether it is operating on DOFs or v-modes. The error, integral, and derivative are all computed in v-mode space:

\[\begin{split}\mathbf{e}^{(v)}_k &= \mathbf{s}[\text{dof\_idx}] - \hat{\mathbf{v}}_k \\ \mathbf{u}^{(v)}_k &= \mathbf{K}_p[\text{dof\_idx}]\,\mathbf{e}^{(v)}_k + \mathbf{K}_i[\text{dof\_idx}]\,\mathbf{I}^{(v)}_k + \mathbf{K}_d[\text{dof\_idx}]\,\mathbf{d}^{(v)}_k\end{split}\]

After the PID produces \(\mathbf{u}^{(v)}_k\), the OFC converts it back to physical DOF corrections:

\[\mathbf{u}_k = -\mathbf{N}\,\mathbf{V}_{h,k}^{T}\,\mathbf{u}^{(v)}_k[\text{:truncate\_index}]\]

This uses get_dofs_from_vmodes(), which takes the first truncate_index elements of the v-mode correction, multiplies by \(\mathbf{V}_{h,k}^T\) to project to normalized DOF space, then de-normalizes with \(\mathbf{N}\).

Important implementation detail: because the PID uses the same gain arrays indexed by dof_idx regardless of whether it operates on DOFs or v-modes, the gain at position \(i\) applies to v-mode \(i\) (ordered by singular value). This means v-mode 0 (highest optical leverage) gets the gain at dof_idx[0], v-mode 1 gets the gain at dof_idx[1], and so on.

Physical meaning of v-mode gains: each entry of \(\mathbf{K}_p\) controls how aggressively a specific optical mode is corrected. For example, assigning a higher gain to v-mode 0 (focus) means the system corrects focus errors more aggressively than higher-order aberrations.

The advantage of v-mode control is the ability to assign different gains based on optical importance:

  • High gain on well-constrained v-modes (large \(\sigma_i\), reliable state estimates) for fast convergence.

  • Low gain on poorly-constrained v-modes (small \(\sigma_i\), noisy estimates) to avoid injecting noise into the corrections.

Equivalent DOF-space gain#

Composing the three steps, the effective DOF-space gain when controlling in v-mode space is:

\[\mathbf{K}_\text{eff} = \mathbf{N}\,\mathbf{V}_h^{T}\,\mathbf{K}^{(v)}_p\, \mathbf{V}_h\,\mathbf{N}^{-1}\]

When the v-mode gain is a uniform scalar (\(\mathbf{K}^{(v)}_p = k_p \mathbf{I}\)):

\[\mathbf{K}_\text{eff} = k_p\,\mathbf{N}\,\mathbf{V}_h^{T}\,\mathbf{V}_h\,\mathbf{N}^{-1} = k_p\,\mathbf{I}\]

The v-mode transformation cancels entirely; v-mode control with uniform gains is identical to DOF-space control with the same scalar gain.

When the v-mode gain has different values per mode, the matrix \(\mathbf{V}_h^{T}\,\mathbf{K}^{(v)}_p\,\mathbf{V}_h\) is a full (non-diagonal) matrix, creating cross-coupling between DOFs. This is the desired behavior: correcting an optical eigenmode (e.g., focus) requires coordinated adjustments to multiple physical DOFs.

For example, with the two-level gain structure \(k_p^{(v)} = 0.18\) for v-modes \(0\)–\(10\) and \(k_p^{(v)} = 0.045\) for v-modes \(11\)–\(21\), the gain decomposes as:

\[\mathbf{K}^{(v)}_p = 0.045\,\mathbf{I} + 0.135\,\mathbf{P}_{11}\]

where \(\mathbf{P}_{11}\) is the projection onto the subspace of the first 11 v-modes. This applies a uniform base gain plus a boost on the most optically important modes.

Correction aggregation#

After each PID step, the DOF correction \(\mathbf{u}_k\) is accumulated into the telescope state:

\[\mathbf{x}_{k+1} = \mathbf{x}_k + \mathbf{u}_k\]

This aggregated state \(\mathbf{x}_{k+1}\) is then decomposed into physical corrections for the individual components:

  • M2 Hexapod: DOFs 0–4 (dZ, dX, dY, rX, rY), transformed via a rotation matrix to hexapod command coordinates.

  • Camera Hexapod: DOFs 5–9, same transformation.

  • M1M3 Bending Modes: DOFs 10–29, converted from bending mode coefficients to actuator forces via the bending-mode-to-force matrix.

  • M2 Bending Modes: DOFs 30–49, same conversion.

These corrections are then issued to the respective subsystem CSCs by ts_mtaos.

Summary#

The following table summarizes the key properties of DOF-space and v-mode control:

Property

DOF space

V-mode space

State vector

Physical DOFs (hexapod positions, bending modes)

Optical eigenmodes ordered by singular value

Gain meaning

Correction fraction per physical DOF

Correction fraction per optical mode

Cross-coupling

None (diagonal gain matrix)

Coordinated multi-DOF corrections

Advantage

Simple, physically intuitive

Gain tuning by optical importance

Number of modes

\(n\) = number of used DOFs

\(n\) (same; from SVD of \(m \times n\) matrix)

Truncation

Applied during state estimation only

Same; PID still operates on all \(n\) modes