Active Suppression of QML

Instead of designing the laser for stable mode-locking, we want to investigate how we could apply state feedback to the pump power to stabilize it.

We start from linearized laser-dynamics:



In [1]:

    
from IPython.display import display
import sympy as sym
sym.init_printing(use_latex='mathjax')
w0, zeta, rho, Toc = sym.symbols('omega_0 zeta rho T_{oc}', real=True)
x = sym.MatrixSymbol('x', 2, 1)
xdot = sym.MatrixSymbol('xdot', 2, 1)
u = sym.MatrixSymbol('u', 1, 1)

A = sym.Matrix([[-2 * w0 * zeta, -w0],
                [w0, 0, ]])
B = sym.Matrix([[rho * w0],
                [0]])
xdot_ = sym.Eq(xdot, A * x + B * u)
xdot_









    Out[1]:





$$\dot{x} = \left[\begin{matrix}\omega_{0} \rho\\0\end{matrix}\right] u + \left[\begin{matrix}- 2 \omega_{0} \zeta & - \omega_{0}\\\omega_{0} & 0\end{matrix}\right] x$$

where

$ \vec{x} = \left[ \begin{array}{ccc} \delta \dot{P} / \omega_0 \\ \delta P \end{array}\right] $ and $ u = \left[ \delta P_P \right] $. $\rho$ is the internal slope efficiency, $\omega_0$ is the natural frequency, and $\zeta$ is the damping ratio of the system.

The output equation is chosen such that it provides the normal output (i.e. $\delta P_{out}$) and additionally the full state:



In [2]:

    
y = sym.MatrixSymbol('y', 3, 1)
C = sym.Matrix([[0, Toc],
                [1, 0],
                [0, 1]])
y_ = sym.Eq(y, C * x)
y_









    Out[2]:





$$y = \left[\begin{matrix}0 & T_{{oc}}\\1 & 0\\0 & 1\end{matrix}\right] x$$

where $T_{oc}$ is the output-coupling ratio.

Pole Placement

(figure from Åström & Murray, "Feedback Systems", figure 6.5)

Now we apply feedback (almost) as shown in the figure. This means that

$ u = k_r \cdot r - K \cdot \vec{y} $

where $r$ is the control signal, $k_r$ is a constant, and $K = \left[0, k_1, k_2 \right]$ is the feedback matrix. Inserting this into the state equations yields:

$ \dot{\vec{x}} = (A - B K C) \vec{x} + B k_r r $

which is a new system $\dot{\vec{x}} = A' \vec{x} + B' r$ with

$ \begin{align*} A' &= A - B K C \\ B' &= B k_r \end{align*} $



In [3]:

    
kr = sym.Symbol('k_r')
r = sym.MatrixSymbol('r', 1, 1)
Ap = sym.MatrixSymbol("A'", 2, 2)
Bp = sym.MatrixSymbol("B'", 2, 1)

k1, k2 = sym.symbols('k_1 k_2', real=True)
K = sym.Matrix([[0, k1, k2]])

Ap_ = sym.Eq(Ap, A - B * K * C)
Bp_ = sym.Eq(Bp, B * kr)

display(Ap_, Bp_)









    





$$A' = \left[\begin{matrix}- k_{1} \omega_{0} \rho - 2 \omega_{0} \zeta & - k_{2} \omega_{0} \rho - \omega_{0}\\\omega_{0} & 0\end{matrix}\right]$$






    





$$B' = \left[\begin{matrix}k_{r} \omega_{0} \rho\\0\end{matrix}\right]$$

Our goal is to to choose $K$ so as to obtain arbitrary poles in the new system. The eigenvalues of $A'$ are:



In [4]:

    
ev = Ap_.rhs.eigenvects()
vals = [v[0] for v in ev]

gamma = sym.Symbol("\gamma", real=True)
gamma_ = sym.Eq(gamma, k1 * rho / 2 + zeta)

nu = sym.Symbol('\\nu', real=True)
nu_ = sym.Eq(nu, (k1 * rho)**2 / 4 + rho * (k1 * zeta - k2) + zeta**2 - 1)

vals = [sym.simplify(v.subs([(k2, sym.solve(nu_, k2)[0]),
                             (zeta, sym.solve(gamma_, zeta)[0])]))
                                                            for v in vals]
display(*vals)









    





$$- \gamma \omega_{0} - \sqrt{\nu} \left\lvert{\omega_{0}}\right\rvert$$






    





$$- \gamma \omega_{0} + \sqrt{\nu} \left\lvert{\omega_{0}}\right\rvert$$

where



In [5]:

    
display(gamma_, nu_)









    





$$\gamma = \frac{k_{1} \rho}{2} + \zeta$$






    





$$\nu = \frac{k_{1}^{2} \rho^{2}}{4} + \rho \left(k_{1} \zeta - k_{2}\right) + \zeta^{2} - 1$$

To obtain the poles as advertised above, $K$ has to be chosen as follows:



In [6]:

    
k1_ = sym.Eq(k1, sym.solve(gamma_, k1)[0])
k2_ = sym.Eq(k2, sym.solve(nu_.subs(k1_.lhs, k1_.rhs), k2)[0])
sym.Eq(sym.MatrixSymbol('K', 1, 3), 
       K.subs(k1_.lhs, k1_.rhs).subs(k2_.lhs, k2_.rhs))









    Out[6]:





$$K = \left[\begin{matrix}0 & \frac{1}{\rho} \left(2 \gamma - 2 \zeta\right) & \frac{1}{\rho} \left(\gamma^{2} - \nu - 1\right)\end{matrix}\right]$$

Finally, we require that the unit-step response on $\delta P$ is $\rho$, i.e. the same as in the uncontrolled system. To achieve that, we choose $k_r$ as follows:



In [7]:

    
dPstep = sym.simplify(-Ap_.rhs**-1 * Bp_.rhs)[1]
dPstep = dPstep.subs(k2_.lhs, k2_.rhs)
kr_ = sym.Eq(kr, sym.solve(sym.Eq(rho, dPstep), kr)[0])
kr_









    Out[7]:





$$k_{r} = \gamma^{2} - \nu$$