clspc

Submodules

Classes

CL_SPC

This class implements the closed-loop subspace predictive control algorithm.

Package Contents

class clspc.CL_SPC(past_window: int, future_window: int, forgetting_factor: float, error_cost: numpy.typing.NDArray[numpy.float64], input_cost: numpy.typing.NDArray[numpy.float64], input_delta_cost: numpy.typing.NDArray[numpy.float64], directional_forgetting: bool = True)[source]

This class implements the closed-loop subspace predictive control algorithm.

Parameters:

  • past_window (int) – Length of past data window.

  • future_window (int) – Length of future data window.

  • forgetting_factor (float) – Forgetting factor.

  • error_cost (npt.NDArray[np.float64]) – Weighting of the reference error term in the cost function.

  • input_cost (npt.NDArray[np.float64]) – Weighting of the input magnitude term in the cost function.

  • input_delta_cost (npt.NDArray[np.float64]) – Weighting of the input change term in the cost function.

  • directional_forgetting (bool) – Enable/disable directional forgetting when calculating the Markov parameters. Defaults to True.

past_window

Length of past data window.

future_window

Length of future data window.

forgetting_factor

Forgetting factor.

q

Weighting of the reference error term.

r

Weighting of the input magnitude term.

r_delta

Weighting of the input change term.

directional_forgetting

Enable/disable directional forgetting in the update of the Markov parameters.

s_delta

\(S_\Delta\). Constant matrix used to calculate input change.

s_v

\(S_v\). Constant matrix used to calculate input change.

phi

\(\varphi_k = \begin{bmatrix} u_{k-p} & y_{k-p} & u_{k-p+1} & y_{k-p+1} & \dots & u_{k-1} & y_{k-1} & \vert & u_{k} & y_{k} \end{bmatrix}^T\)

Note

self.phi includes \(y_{k}\) for convenience in the code. \(y_{k}\) is not included in \(\varphi_k\) for calculations.

markov_parameters

\(\hat{\Theta} = \begin{bmatrix} \Xi^{(u_{k-p})} & \Xi^{(y_{k-p})} & \Xi^{(u_{k-p+1})} & \Xi^{(y_{k-p+1})} & \dots & \Xi^{(u_{k-1})} & \Xi^{(y_{k-1})} & \vert & \Xi^{(u_{k})} \end{bmatrix}^T\)

cholesky_factor

Cholesky factorisation of the covariance matrix.

steps = 0

Step counter.

calculate_s_delta(f: int) → numpy.typing.NDArray[numpy.float64][source]

Calculate constant matrix \(S_\Delta\).

\[\begin{split}S_\Delta \in \mathbb{R}^{f\times f} = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ -1 & 1 & 0 & \dots & 0 \\ 0 & -1 & 1 & & \vdots \\ \vdots & & \ddots & \ddots & 0 \\ 0 & \dots & 0 & -1 & 1 \end{bmatrix}\end{split}\]
Parameters:

f (int) – Length of future window.

Returns:

Matrix \(S_\Delta\).

Return type:

npt.NDArray[np.float64]

calculate_s_v(p: int, f: int) → numpy.typing.NDArray[numpy.float64][source]

Calculate constant matrix \(S_v\).

\[\begin{split}S_v \in \mathbb{R}^{f\times p} = \begin{bmatrix} 0 & \dots & 0 & 1 \\ 0 & \dots & 0 & 0 \\ \vdots & & \vdots & \vdots \\ 0 & \dots & 0 & 0 \end{bmatrix}\end{split}\]
Parameters:

  • p (int) – Length of past window.

  • f (int) – Length of future window.

Returns:

Matrix \(S_v\).

Return type:

npt.NDArray[np.float64]

step(uk: float, yk: float, rf: numpy.typing.NDArray[numpy.float64]) → float[source]

Perform one control step:

  1. Update self.phi with \(u_k\), \(y_k\)

  2. Increase step counter

If the step counter is bigger than the past window size:

  1. Update Markov parameters

  2. Construct \(\tilde{\mathcal{H}}\)

  3. Calculate \(\mathcal{G}\), \(\Gamma \tilde{\mathcal{K}}\), \(\Gamma \tilde{\mathcal{L}}\)

  4. Solve quadratic program

  5. Return \(u_{k+1}\).

Parameters:

  • uk (float) – \(u_k\), most recent input.

  • yk (float) – \(y_k\), most recent output.

  • rf (npt.NDArray[np.float64]) – \(r_f\), future reference array of size CL_SPC.future_window.

Returns:

\(u_{k+1}\), new input.

Return type:

float

update_markov_parameters(phi: numpy.typing.NDArray[numpy.float64], yk: float) → numpy.typing.NDArray[numpy.float64][source]

Update Markov parameters:

  1. Create pre-array

  2. Perform directional forgetting step using Givens rotations

  3. Perform covariance update using Givens rotations

  4. Update Markov parameters

Parameters:

  • phi (npt.NDArray[np.float64]) – self.phi, not including \(y_k\).

  • yk (float) – System output \(y_k\).

Returns:

Updated Markov parameters \(\hat{\Theta}_k\)

Return type:

npt.NDArray[np.float64]

static perform_givens_rotations(pre_array: numpy.typing.NDArray[numpy.float64], row_range: Iterable[int], r_col_range: Iterable[int], zero_col_range: Iterable[int]) → numpy.typing.NDArray[numpy.float64][source]

Perform a set of givens rotations to zero each element bi = pre_array[row_range[i], zero_col_range[i]] by rotating it to onto the corresponding ai = pre_array[row_range[i], r_col_range[i]].

\[\begin{split}\begin{bmatrix} a & b \end{bmatrix} \begin{bmatrix} c & -s \\ s & c \end{bmatrix} = \begin{bmatrix} r & 0 \end{bmatrix}\end{split}\]
\[r = \sqrt{a^2 + b^2}\]
\[c = \frac{a}{r}\]
\[s = \frac{b}{r}\]
Parameters:

  • pre_array (npt.NDArray[np.float64]) – Input array containing elements to be zeroed.

  • row_range (Iterable[int]) – _description_

  • r_col_range (Iterable[int]) – _description_

  • zero_col_range (Iterable[int]) – _description_

Returns:

_description_

Return type:

npt.NDArray[np.float64]

calculate_h_tilde(markov_y: numpy.typing.NDArray[numpy.float64]) → numpy.typing.NDArray[numpy.float64][source]

Construct \(\tilde{\mathcal{H}}\) from Markov parameters.

\[\begin{split}\tilde{\mathcal{H}} = \begin{bmatrix} I & 0 & \dots & 0 \\ -CK & I & & \vdots \\ \vdots & \vdots & \ddots & 0 \\ - C\tilde{A}^{f-2}K & -C\tilde{A}^{f-3}K & \dots & I \end{bmatrix}\end{split}\]

with

\[\tilde{\Xi}^{(y_{k-i})} = C\tilde{A}^{i-1}K\]
Parameters:

markov_y (npt.NDArray[np.float64]) – \(\tilde{\Xi}^{(y)}\), Markov parameters pertaining to the outputs.

Returns:

\(\tilde{\mathcal{H}}\).

Return type:

npt.NDArray[np.float64]

calculate_g(markov_u: numpy.typing.NDArray[numpy.float64], h_tilde: numpy.typing.NDArray[numpy.float64]) → numpy.typing.NDArray[numpy.float64][source]

Calculate \(\mathcal{G}\). First, construct \(\tilde{\mathcal{G}}\) from Markov parameters, then calculate \(\mathcal{G}\) using \(\tilde{\mathcal{H}}\) taking advantage of the fact that the latter is lower triangular.

\[\begin{split}\tilde{\mathcal{G}} = \begin{bmatrix} D & 0 & \dots & 0 \\ C\tilde{B} & D & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ C\tilde{A}^{f-2}\tilde{B} & C\tilde{A}^{f-3}\tilde{B} & \dots & D \end{bmatrix}\end{split}\]

with

\[\begin{split}\tilde{\Xi}^{(u_{k-i})} = \begin{cases} D, &\text{if } i = 0 \\ C\tilde{A}^{i-1}\tilde{B}, &\text{if } i > 0 \end{cases}\end{split}\]

Then

\[\mathcal{G} = \tilde{\mathcal{H}}^{-1} \tilde{\mathcal{G}}\]
Parameters:

  • markov_u (npt.NDArray[np.float64]) – \(\tilde{\Xi}^{(u)}\), Markov parameters pertaining to the inputs.

  • h_tilde (npt.NDArray[np.float64]) – \(\tilde{\mathcal{H}}\).

Returns:

\(\mathcal{G}\).

Return type:

npt.NDArray[np.float64]

calculate_gamma_l(markov_u: numpy.typing.NDArray[numpy.float64], h_tilde: numpy.typing.NDArray[numpy.float64]) → numpy.typing.NDArray[numpy.float64][source]

Calculate \(\Gamma \tilde{\mathcal{L}}\). First, construct \(\tilde{\Gamma}\tilde{\mathcal{L}}\) from Markov parameters, then calculate \(\Gamma \tilde{\mathcal{L}}\) using \(\tilde{\mathcal{H}}\) taking advantage of the fact that the latter is lower triangular.

\[\begin{split}\tilde{\Gamma}\tilde{\mathcal{L}} = \begin{bmatrix} \tilde{\Xi}^{(u_{k-p})} & \tilde{\Xi}^{(u_{k-p+1})} & \dots & \tilde{\Xi}^{(u_{k-p+f-1})} & \dots & \tilde{\Xi}^{(u_{k-1})} \\ 0 & \tilde{\Xi}^{(u_{k-p})} & \dots & \tilde{\Xi}^{(u_{k-p+f-2})} & \dots & \tilde{\Xi}^{(u_{k-2})} \\ \vdots & & & \vdots & \ddots & \vdots \\ 0 & & \dots & \tilde{\Xi}^{(u_{k-p})} & \dots & \tilde{\Xi}^{(u_{k-f})} \end{bmatrix}\end{split}\]

Then

\[\Gamma \tilde{\mathcal{L}} = \tilde{\mathcal{H}}^{-1} \tilde{\Gamma} \tilde{\mathcal{L}}\]
Parameters:

  • markov_u (npt.NDArray[np.float64]) – \(\tilde{\Xi}^{(u)}\), Markov parameters pertaining to the inputs.

  • h_tilde (npt.NDArray[np.float64]) – \(\tilde{\mathcal{H}}\).

Returns:

\(\Gamma \tilde{\mathcal{L}}\).

Return type:

npt.NDArray[np.float64]

calculate_gamma_k(markov_y: numpy.typing.NDArray[numpy.float64], h_tilde: numpy.typing.NDArray[numpy.float64]) → numpy.typing.NDArray[numpy.float64][source]

Calculate \(\Gamma \tilde{\mathcal{K}}\). First, construct \(\tilde{\Gamma}\tilde{\mathcal{K}}\) from Markov parameters, then calculate \(\Gamma \tilde{\mathcal{K}}\) using \(\tilde{\mathcal{H}}\) taking advantage of the fact that the latter is lower triangular.

\[\begin{split}\tilde{\Gamma}\tilde{\mathcal{K}} = \begin{bmatrix} \tilde{\Xi}^{(y_{k-p})} & \tilde{\Xi}^{(y_{k-p+1})} & \dots & \tilde{\Xi}^{(y_{k-p+f-1})} & \dots & \tilde{\Xi}^{(y_{k-1})} \\ 0 & \tilde{\Xi}^{(y_{k-p})} & \dots & \tilde{\Xi}^{(y_{k-p+f-2})} & \dots & \tilde{\Xi}^{(y_{k-2})} \\ \vdots & & & \vdots & \ddots & \vdots \\ 0 & & \dots & \tilde{\Xi}^{(y_{k-p})} & \dots & \tilde{\Xi}^{(y_{k-f})} \end{bmatrix}\end{split}\]

Then

\[\Gamma \tilde{\mathcal{K}} = \tilde{\mathcal{H}}^{-1} \tilde{\Gamma} \tilde{\mathcal{K}}\]
Parameters:

  • markov_y (npt.NDArray[np.float64]) – \(\tilde{\Xi}^{(y)}\), Markov parameters pertaining to the outputs.

  • h_tilde (npt.NDArray[np.float64]) – \(\tilde{\mathcal{H}}\).

Returns:

\(\Gamma \tilde{\mathcal{K}}\).

Return type:

npt.NDArray[np.float64]

_calculate_gamma_kl(markov: numpy.typing.NDArray[numpy.float64], h_tilde: numpy.typing.NDArray[numpy.float64]) → numpy.typing.NDArray[numpy.float64][source]

Implements common functionality of CL_SPC.calculate_gamma_k() and CL_SPC.calculate_gamma_l(), since array structure is identical and only differs in which Markov parameters are used.

Parameters:

  • markov (npt.NDArray[np.float64]) – \(\tilde{\Xi}^{(u)}\) or \(\tilde{\Xi}^{(y)}\), Markov parameters.

  • h_tilde (npt.NDArray[np.float64]) – \(\tilde{\mathcal{H}}\).

Returns:

\(\Gamma \tilde{\mathcal{K}}\) or \(\Gamma \tilde{\mathcal{L}}\), depending on provided Markov parameters.

Return type:

npt.NDArray[np.float64]

solve_qp(up: numpy.typing.NDArray[numpy.float64], yp: numpy.typing.NDArray[numpy.float64], rf: numpy.typing.NDArray[numpy.float64], g: numpy.typing.NDArray[numpy.float64], gamma_l: numpy.typing.NDArray[numpy.float64], gamma_k: numpy.typing.NDArray[numpy.float64]) → float[source]

Solve quadratic program

\[J = u^T_f P u_f + q^T u_f\]

with

\[P = \mathcal{G}^T Q_a \mathcal{G} + S^T_\Delta R^\Delta_a S_\Delta + R_a\]
\[q^T = 2\left(u^T_p \Gamma \tilde{\mathcal{L}} Q_a \mathcal{G} + u^T_p \Gamma \tilde{\mathcal{K}} Q_a \mathcal{G} - r^T_f Q \mathcal{G} - u^T_p S^T_v R^\Delta_a S_\Delta \right)\]
Parameters:

  • up (npt.NDArray[np.float64]) – \(u_p\), array of past inputs of length CL_SPC.past_window.

  • yp (npt.NDArray[np.float64]) – \(y_p\), array of past outputs of length CL_SPC.past_window.

  • rf (npt.NDArray[np.float64]) – \(r_f\), array of future references of length CL_SPC.future_window.

  • g (npt.NDArray[np.float64]) – \(\mathcal{G}\).

  • gamma_l (npt.NDArray[np.float64]) – \(\Gamma \tilde{\mathcal{L}}\).

  • gamma_k (npt.NDArray[np.float64]) – \(\Gamma \tilde{\mathcal{K}}\).

Returns:

\(u_{k+1}\), next input.

Return type:

float