Skip to content

qttools.boundary_conditions.lyapunov.doubling#

[docs] module qttools.boundary_conditions.lyapunov.doubling 

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
# Copyright (c) 2024-2026 ETH Zurich and the authors of the qttools package.

import warnings

from qttools import NDArray, xp
from qttools.boundary_conditions.lyapunov.lyapunov import LyapunovSolver


class Doubling(LyapunovSolver):
    """A solver for the Lyapunov equation using iterative doubling.

    Parameters
    ----------
    max_iterations : int, optional
        The maximum number of iterations to perform.
    convergence_rel_tol : float, optional
        The required relative accuracy for convergence.
    convergence_abs_tol : float, optional
        The required absolute accuracy for convergence.
        Either convergence_rel_tol or convergence_abs_tol must be satisfied.

    """

    def __init__(
        self,
        max_iterations: int = 10,
        convergence_rel_tol: float = 1e-5,
        convergence_abs_tol: float = 1e-8,
    ) -> None:
        """Initializes the solver."""
        self.max_iterations = max_iterations
        self.convergence_rel_tol = convergence_rel_tol
        self.convergence_abs_tol = convergence_abs_tol

    def __call__(
        self,
        a: NDArray,
        q: NDArray,
        contact: str,
    ) -> NDArray:
        """Computes the solution of the discrete-time Lyapunov equation.

        The matrices a and q can have different ndims with q.ndim >= a.ndim (will broadcast)

        Parameters
        ----------
        a : NDArray
            The system matrix.
        q : NDArray
            The right-hand side matrix.
        contact : str
            The contact to which the boundary blocks belong.

        Returns
        -------
        x : NDArray
            The solution of the discrete-time Lyapunov equation.

        """

        assert q.shape[-2:] == a.shape[-2:]
        assert q.ndim >= a.ndim

        a = xp.broadcast_to(a, q.shape)
        a_i = a.copy()
        x = q.copy()

        for __ in range(self.max_iterations):
            x_i = x + a_i @ x @ a_i.conj().swapaxes(-1, -2)

            absolute_recursion_errors = xp.linalg.norm(x_i - x, axis=(-2, -1))
            relative_recursion_errors = absolute_recursion_errors / xp.linalg.norm(
                x_i, axis=(-2, -1)
            )
            x = x_i

            if xp.all(
                (relative_recursion_errors < self.convergence_rel_tol)
                | (absolute_recursion_errors < self.convergence_abs_tol)
            ):
                break

            a_i = a_i @ a_i

        else:  # Did not break, i.e. max_iterations reached.
            warnings.warn(
                "The doubling method for the Lyapunov equation did not converge.",
                RuntimeWarning,
            )

        return x