diff --git a/content/zettels/feedforward_control.md b/content/zettels/feedforward_control.md
index dc25277..702afce 100644
--- a/content/zettels/feedforward_control.md
+++ b/content/zettels/feedforward_control.md
@@ -105,14 +105,14 @@ and \\(s\\) the snap, \\(j\\) the jerk, \\(a\\) the acceleration and \\(v\\) the
 
 The same architecture shown in Figure [3](#figure--fig:feedforward-fourth-order-feedforward-architecture) can be used.
 
-In order to implement a fourth order trajectory, look at [this](https://www.mathworks.com/matlabcentral/fileexchange/16352-advanced-setpoints-for-motion-systems) nice implementation in Simulink of fourth-order trajectory planning (see also (<a href="#citeproc_bib_item_1">Lambrechts, Boerlage, and Steinbuch 2004</a>)).
+In order to implement a fourth order trajectory, look at [this](https://www.mathworks.com/matlabcentral/fileexchange/16352-advanced-setpoints-for-motion-systems) nice implementation in Simulink of fourth-order trajectory planning (see also <&lambrechts04_trajec>).
 
 
 ## Model Based Feedforward Control for Second Order resonance plant {#model-based-feedforward-control-for-second-order-resonance-plant}
 
 <span class="org-target" id="org-target--sec-model-based-feedforward"></span>
 
-See (<a href="#citeproc_bib_item_2">Schmidt, Schitter, and Rankers 2020</a>) (Section 4.2.1).
+See <&schmidt20_desig_high_perfor_mechat_third_revis_edition> (Section 4.2.1).
 
 Suppose we have a second order plant (could typically be a piezoelectric stage):
 \\[ G(s) = \frac{C\_f \omega\_0^2}{s^2 + 2\xi \omega\_0 s + \omega\_0^2} \\]
@@ -136,7 +136,90 @@ For \\(\xi = 1\\) this results in the following transfer function:
 
 {{< figure src="/ox-hugo/feedforward_compensated_system.png" caption="<span class=\"figure-number\">Figure 5: </span>Bode plot of the feedforward controlled system" >}}
 
-<style>.csl-entry{text-indent: -1.5em; margin-left: 1.5em;}</style><div class="csl-bib-body">
-  <div class="csl-entry"><a id="citeproc_bib_item_1"></a>Lambrechts, P., M. Boerlage, and M. Steinbuch. 2004. “Trajectory Planning and Feedforward Design for High Performance Motion Systems.” In <i>Proceedings of the 2004 American Control Conference</i>, nil. doi:<a href="https://doi.org/10.23919/acc.2004.1384042">10.23919/acc.2004.1384042</a>.</div>
-  <div class="csl-entry"><a id="citeproc_bib_item_2"></a>Schmidt, R Munnig, Georg Schitter, and Adrian Rankers. 2020. <i>The Design of High Performance Mechatronics - Third Revised Edition</i>. Ios Press.</div>
-</div>
+
+## Advanced Feedforward (from MIMO training) {#advanced-feedforward--from-mimo-training}
+
+A typical control configuration for motion systems consists of:
+
+-   A setpoint generator (SPG)
+-   A feedback controller (\\(K\_{fb}\\))
+-   A feedforward controller (\\(K\_{ff}\\))
+
+{{< figure src="/ox-hugo/feedforward_schematic.png" >}}
+
+The closed-loop error (no disturbances) is:
+\\[ e(s) = (1 + G(s)K\_{fb})^{-1} (1 - G(s)K\_{ff}(s)) r(s) \\]
+It therefore depends on:
+
+1.  the setpoint \\(r\\)
+2.  the feedforward controller
+3.  the feedback controller
+
+**Setpoint generation**:
+
+-   It can be 2nd order, 3rd order or 4th order
+-   For 4th order, derivative of jerk is generated over time, and then integrated 4 times to give: jerk, acceleration, velocity and position.
+
+**2nd order setpoint generation**:
+If we compute the fourier transform of the generated acceleration, we get the following signal (-20db/dec):
+![](figs/feedforward_2nd_order_fourier.png)
+Notches are at \\(f\_1\\), \\(2f\_1\\), \\(3f\_1\\), ... with \\(f\_1 = \frac{a\_{\text{max}}}{v\_{\text{max}}}\\).
+It is therefore possible to choose the velocity and acceleration such that \\(f\_1\\) (or one of its integral multiple) matches the resonance frequency of the system.
+Therefore, the acceleration time constant can be chosen at the inverse of the plant resonance.
+
+**3rd order setpoint generation**:
+There is a drawback of having an extra time of \\(\frac{a\_{max}}{J\_{max}}\\) seconds.
+However, we get an additional -20db/dec at high frequency, and additional notches at \\(f\_2 = \frac{j\_{max}}{a\_{max}}\\).
+This new notch has larger "damping" and can be used to be more robust against resonances of the plant.
+
+**Feedforward control**:
+Plant inversion: if \\(K\_{ff} = G^{-1}(s) \Longrightarrow e(s) = 0\\)
+Challenges:
+
+-   Model required
+-   High order
+-   Delay/non-minimum phase?
+
+**Rigid body dynamics**:
+\\(G(s) = \frac{1}{ms^2}\\)
+In that case, \\(G^{-1}(s) = ms^2\\), and with 2nd order setpoint, a feedforward controller \\(K\_{ff}(s) = m\\) gives good performances.
+
+{{< figure src="/ox-hugo/feedforward_schematic_rigid_body.png" >}}
+
+**Discrete time implementation**.
+The DAC can usually be modelled by a "Zero Order Hold" (ZOH) and the ADC with a "sampler".
+This adds **1.5 samples of delay**: \\(0.5z^{-1} + 0.5z^{-2}\\).
+
+It seems the ZOH can be modelled by an "half-sample delay".
+There is an additional one sample delay.
+
+{{< figure src="/ox-hugo/feedforward_schematic_zoh_sampler.png" >}}
+
+Therefore, it is very important to match the delay of the plant:
+
+> In high-performance control system, it can be useful to consider propagation delay when designing the feedforward and feedback controllers.
+> Feedforward control directly uses the reference trajectory and does not depend on any measurement data.
+> However, the feedback controller uses (delayed) measured position data.
+> Due to propagation delay in the control system (caused by the controller, actuator or sensor), it can take multiple cycles for the effect of feedforward control to be observed in the measured position.
+> In the meantime, the feedback control is already seeing a tracking error and is compensating for it.
+> Essentially, the result of the feedforward action arrives too late, resulting in possible overcompensation by the feedback control.
+> When the propagation delay in the control system is known, it can be compensated for by applying this same delay to the demand position in the tracking error calculation.
+
+{{< figure src="/ox-hugo/feedforward_schematic_delay.png" >}}
+
+**Feedforward for flexible dynamics**:
+4th order dynamics:
+
+{{< figure src="/ox-hugo/feedforward_4th_order.png" >}}
+
+Dynamics from \\(F\\) to \\(x\_2\\) is:
+\\[ G(s) = \frac{x\_2}{F} = \frac{cs + k}{m\_1m\_2 s^2(s^2 + 2\xi\omega\_0 s + \omega\_0^2)} \\]
+We take the inverse for the feedforward controller:
+\\[ K\_{ff}(s) = G^{-1}(s) = \frac{m\_1m\_2 s^2(s^2 + 2\xi\omega\_0 s + \omega\_0^2)}{cs + k} \\]
+If we neglect damped: \\(\xi = 0\\), and we get:
+\\[ K\_{ff}(s) = \underbrace{\frac{m\_1m\_2}{k} s^4}\_{\text{snap FF}} + \underbrace{(m\_1 + m\_2) s^2}\_{\text{acc FF}} \\]
+This can be solved by using **snap feedforward**
+
+{{< figure src="/ox-hugo/feedforward_schematic_snap.png" >}}
+
+<./biblio/references.bib>
diff --git a/content/zettels/multivariable_control.md b/content/zettels/multivariable_control.md
index e8e49af..9525203 100644
--- a/content/zettels/multivariable_control.md
+++ b/content/zettels/multivariable_control.md
@@ -10,6 +10,31 @@ Tags
 A very nice book about Multivariable Control is <skogestad07_multiv_feedb_contr>
 
 
+## Transfer functions for Multi-Input Multi-Output systems {#transfer-functions-for-multi-input-multi-output-systems}
+
+{{< figure src="/ox-hugo/mimo_tf.png" >}}
+
+\\[ T\_i = -\frac{u}{d\_i} = (I + KG)^{-1} KG \\]
+\\[ T\_o = -\frac{p\_o}{d\_o} = (I + GK)^{-1} GK \\]
+\\[ S\_i = \frac{p\_i}{d\_i} = (I + KG)^{-1} \\]
+\\[ S\_o = \frac{y}{d\_o} = (I + GK)^{-1} \\]
+
+
+## Measures of interaction {#measures-of-interaction}
+
+-   Interaction index (for \\(2 \times 2\\) plant):
+    \\[ \phi = \frac{g\_{12}g\_{21}}{g\_{11}g\_{22}} \\]
+    When \\(\phi\\) is close to zero, this means there is no interaction.
+-   The **relative gain array** of a square matrix:
+    \\[ \text{RGA}(G) \triangleq G \times ( G^{-1})^T \\]
+
+
+## Stability {#stability}
+
+-   **Characteristic Loci**: Eigenvalues of \\(G(j\omega)\\) plotted in the complex plane
+-   **Generalized Nyquist Criterion**: If \\(G(s)\\) has \\(p\_0\\) unstable poles, then the closed-loop system with return ratio \\(kG(s)\\) is stable if and only if the characteristic loci of \\(kG(s)\\), taken together, encircle the point \\(-1\\), \\(p\_0\\) times anti-clockwise, assuming there are no hidden modes
+
+
 ## Bibliography {#bibliography}
 
 <./biblio/references.bib>
diff --git a/content/zettels/nyquist_stability_criterion.md b/content/zettels/nyquist_stability_criterion.md
new file mode 100644
index 0000000..2626e3d
--- /dev/null
+++ b/content/zettels/nyquist_stability_criterion.md
@@ -0,0 +1,72 @@
++++
+title = "Nyquist stability criterion"
+author = ["Dehaeze Thomas"]
+draft = false
++++
+
+Tags
+:
+
+
+## Theory {#theory}
+
+The main reason why the Nyquist plot is used it that it can be used with the experimental FRF data!
+
+The zeros and pole of a MIMO system are the zeros and pole of the determinant of \\(G(s)\\).
+\\[ \det(G(s)) = \frac{z\_G(s)}{p\_G(s)} \\]
+The polynomial \\(p\_G(s)\\) is normally called the **open-loop characteristic** polynomial.
+
+In a MIMO feedback system:
+
+-   The transfer function matrix open-loop is:
+    \\[ L(s) = G(s) K(s) \neq K(s) G(s) \\]
+-   The transfer function matrix closed-loop is:
+    \\[ T(s) = [I + L(s)]^{-1} L(s) \\]
+-   **Return difference matrix**:
+    \\[ F(s) = [I + L(s)] \\]
+
+The closed-loop system is stable if the zeros of the closed-loop characteristic polynomial lie in the complex open left half plane.
+There are the zeros of:
+\\[ \det(I + GK) \\]
+
+<div class="important">
+
+**MIMO Nyquist stability criteria**:
+\\[ \det(I + G(s)K(s)) = 0 \quad \text{for} \quad \text{Re}(s)<0 \\]
+To check the closed-loop stability graphically, plot the Nyquist of \\(\det(I + GK)\\) and evaluate the encirclement with respect to the point \\((0,0)\\).
+The Nyquist plot is the image of the imaginary axis (\\(j\omega\\)) under \\(\det(I + GK)\\), i.e. it is the evolution of \\(\det(I + G(j\omega)K(j\omega))\\) in the complex plane.
+Note that there is a single plot, even in the MIMO case.
+
+</div>
+
+<div class="important">
+
+**Eigenvalue loci**:
+The eigenvalue loci (sometimes called the characteristic loci) are defined as the eigenvalues of the frequency response function of the open-loop transfer function matrix \\(G(s)K(s)\\).
+This time, there are \\(n\\) plots, where \\(n\\) is the size of the system.
+
+</div>
+
+
+## Matlab Example {#matlab-example}
+
+Sure we have identified a system with 6 inputs and 6 outputs.
+The Matlab object has dimension `6 x 6 x n` with `n` is the number of frequency points.
+
+First, compute the open-loop gain:
+
+```matlab
+L  = zeros(6, 6, length(f));
+
+for i_f = 1:length(f)
+    L(:,:,i_f)  = squeeze(G(:,:,i_f))*freqresp(K, f(i_f), 'Hz');
+end
+```
+
+Then, compute the eigenvalues of this open-loop gain:
+Finally, plot the (complex) eigenvalues in the complex plane:
+
+
+## Bibliography {#bibliography}
+
+<./biblio/references.bib>
diff --git a/content/zettels/relative_gain_array.md b/content/zettels/relative_gain_array.md
new file mode 100644
index 0000000..232ab62
--- /dev/null
+++ b/content/zettels/relative_gain_array.md
@@ -0,0 +1,41 @@
++++
+title = "Relative Gain Array"
+author = ["Dehaeze Thomas"]
+draft = false
++++
+
+Tags
+:
+
+Imagine a 2x2 plant:
+\\[ G(s) = \begin{bmatrix} g\_{11}(s) & g\_{12}(s)\\\ g\_{21}(s) & g\_{22}(s) \end{bmatrix} \\]
+
+Now suppose a controller \\(k\_{1}(s)\\) is closed on the first loop.
+The new transfer function from \\(u\_2\\) to \\(y\_2\\) (while the first loop is closed) is :
+
+\begin{equation}
+y\_2 = g\_{22} \left[ 1 - \underbrace{\frac{g\_{21}g\_{12}}{g\_{11}g\_{22}}}\_{\phi} \underbrace{\frac{g\_{11} k\_{11}}{1  + g\_{11}k\_{11}}}\_{T\_{11}} \right]
+\end{equation}
+
+\\(\phi\\) is called the **interaction index**.
+\\(T\_{11}\\) is the complementary sensitivity of the first loop (equal to 1 in the bandwidth of the first controller).
+Therefore, we want \\(\phi \approx 0\\) to have no interaction.
+
+Similarly, the **relative gain** is defined as:
+\\[ \Lambda = \frac{1}{1 - \phi} \\]
+And we want \\(\Lambda \approx 1\\) to have no interaction.
+
+Note that the scaling of inputs or outputs of the MIMO plant has no effect on \\(\phi\\) or \\(\Lambda\\).
+
+The **relative gain array** is defined as:
+\\[ \Lambda(G) = G \star (G^{-1})^{T} \\]
+where \\(\star\\) means element wise multiplication.
+
+```matlab
+RGA = G.*pinv(G.')
+```
+
+
+## Bibliography {#bibliography}
+
+<./biblio/references.bib>
diff --git a/static/ox-hugo/feedforward_4th_order.png b/static/ox-hugo/feedforward_4th_order.png
new file mode 100644
index 0000000..dc343d5
Binary files /dev/null and b/static/ox-hugo/feedforward_4th_order.png differ
diff --git a/static/ox-hugo/feedforward_schematic.png b/static/ox-hugo/feedforward_schematic.png
new file mode 100644
index 0000000..b73dc1c
Binary files /dev/null and b/static/ox-hugo/feedforward_schematic.png differ
diff --git a/static/ox-hugo/feedforward_schematic_delay.png b/static/ox-hugo/feedforward_schematic_delay.png
new file mode 100644
index 0000000..90c4d09
Binary files /dev/null and b/static/ox-hugo/feedforward_schematic_delay.png differ
diff --git a/static/ox-hugo/feedforward_schematic_rigid_body.png b/static/ox-hugo/feedforward_schematic_rigid_body.png
new file mode 100644
index 0000000..5d662ba
Binary files /dev/null and b/static/ox-hugo/feedforward_schematic_rigid_body.png differ
diff --git a/static/ox-hugo/feedforward_schematic_snap.png b/static/ox-hugo/feedforward_schematic_snap.png
new file mode 100644
index 0000000..3b96f7e
Binary files /dev/null and b/static/ox-hugo/feedforward_schematic_snap.png differ
diff --git a/static/ox-hugo/feedforward_schematic_zoh_sampler.png b/static/ox-hugo/feedforward_schematic_zoh_sampler.png
new file mode 100644
index 0000000..2c014d0
Binary files /dev/null and b/static/ox-hugo/feedforward_schematic_zoh_sampler.png differ
diff --git a/static/ox-hugo/mimo_tf.png b/static/ox-hugo/mimo_tf.png
new file mode 100644
index 0000000..ac205ee
Binary files /dev/null and b/static/ox-hugo/mimo_tf.png differ