Notebook

Automatic Matching¶

This notebooks explains how the automatic matching is implemented on WEST ICRH antenna.

In [ ]:

%load_ext autoreload
%autoreload 2

In [ ]:

import matplotlib.pyplot as plt
import skrf as rf
import numpy as np
import sys
from tqdm.notebook import tqdm

sys.path.append('..')
from west_ic_antenna import WestIcrhAntenna

Introduction¶

The WEST automatic matching of the capacitors of an ICRH antenna is made with a negative feedback loop. ICRH Operator defines the desired input impedance at the T-junction setpoint $Z_{T,SP}$. This setpoint is compared to the actual impedance at the T-junction $Z_T$ and capacitors are adjusted to minimize this difference in realtime. The input impedance at the T-junction is determined from the input impedance at the bidirective coupler located behind the antennas.

The different steps are described below:

RF Measurements: the forward $P_{fwd}$ and reflected powers $P_{ref}$ are measured at the bidirective coupler located behind an antenna. The phase difference between these two $\Phi_\Gamma$ is also measured. The impedance at the bidirective coupler $Z_{coupler}$ is deduced from the complex reflection coefficient $\Gamma$ :

$$ \Gamma = \sqrt{\frac{P_{ref}}{P_{fwd}}} \exp⁡ \left[ j(\Phi_\Gamma + \Phi_c ) \right]$$$$ Z_{coupler} = Z_0 \frac{1 + \Gamma}{1 - \Gamma} $$

where $Z_0$=30 Ω is the characteristic impedance of the transmission line and $\Phi_c$ a possible phase correction that the IC operator can tune.

$Z_T$ Calculation:

Let $A$ the ABCD matrix of the assembly of the piece of transmission line between the bidirective coupler and the antenna elements up to the T-junction. This matrix has been deduced from the CAD model of the antenna ($S_{SWIT}$) and from the measurement of the piece of transmission lines ($S_{ciseau}$) for each antennas (which can differ in length).

The impedance at the T-junction is deduced from the impedance $Z_{coupler}$ (here asssuming port 1 is located at the T junction):

$$ Z_T = \frac{A_{11} Z_{coupler} + A_{12}}{A_{21} Z_{coupler} + A_{22} } $$

NB: if port 1 is instead the 30 Ohm side, the equation would have been:

$$ Z_T = \frac{A_{12} - Z_{coupler} A_{22}}{A_{21} Z_{coupler} - A_{11} } $$

Normalized error signal: the relative change of the measured error $\varepsilon$ is :

$$ \mathbf{\varepsilon} = \left( \begin{array}{c} \Re(\varepsilon) \\ \Im(\varepsilon) \end{array} \right) = \left( \begin{array}{c} \Re\left( \frac{Z_{T,sp}}{Z_T} \right) - 1 \\ \Im\left( \frac{Z_{T,sp}}{Z_T} \right) \end{array} \right) $$

Controller: gives a velocity request $V_r=(V_{r,top}, V_{r,bot})$ to the capacitors actuator:

$$ \left( \begin{array}{c} V_{r,top} \\ V_{r,bot} \end{array} \right) = \mathbf{T} \times \mathbf{\varepsilon} $$

The derivation of the matrix $\mathbf{T}$ is detailed in the next section

Controler Development¶

Introduction¶

The impedance at the T-junction $Z_T$ is a function of the capacitances set $\mathbf{C}=( C_{top}, C_{bot})^t$:

$$ Z_T = Z_T(\mathbf{C}) $$

Close to the desired capacitances solution $\mathbf{C}_{SP}=(C_{top,SP}, C_{bot,SP})^t$, the previous equation can be approximated thanks to the Taylor theorem:

$$ Z_T(\mathbf{C}) \approx Z_T(\mathbf{C}_{SP}) + \nabla Z_T(\mathbf{C}_{SP}) \cdot (\mathbf{C} - \mathbf{C}_{SP}) $$

with $$ \nabla Z_T(\mathbf{C}_{SP}) = \left. \left( \begin{array}{c} \frac{\partial Z_T}{\partial C_{top}} \\ \frac{\partial Z_T}{\partial C_{bot}} \end{array} \right) \right|_{\mathbf{C}=\mathbf{C}_{SP}} $$ We define later $\delta\mathbf{C}=\mathbf{C} - \mathbf{C}_{SP}$ and recalls that the first term on the right hand side $Z_T(\mathbf{C}_{SP})$ corresponds to the desired Set Point $Z_{T,SP}$.

Splitting $Z_T$ into real and imaginary parts leads to:

$$ \left( \begin{array}{c} \Re(Z_T)(\mathbf{C}) \\ \Im(Z_T)(\mathbf{C}) \end{array} \right) \approx \left( \begin{array}{c} \Re(Z_{T,SP})\\ \Im(Z_{T,SP}) \end{array} \right) + \mathbf{D} \cdot \delta\mathbf{C} $$

where $$ \mathbf{D} = \left. \left( \begin{array}{cc} \frac{\partial \Re(Z_T)}{\partial C_{top}} & \frac{\partial \Re(Z_T)}{\partial C_{bot}} \\ \frac{\partial \Im(Z_T)}{\partial C_{top}} & \frac{\partial \Im(Z_T)}{\partial C_{bot}} \end{array} \right) \right|_{\mathbf{C}=\mathbf{C}_{opt}} $$

The previous formula can be inversed to given the required $\delta \mathbf{C}$ increment to reach the Set Point:

$$ \delta\mathbf{C} \approx \mathbf{D}^{-1} \left( \left( \begin{array}{c} \Re(Z_T)(\mathbf{C}) \\ \Im(Z_T)(\mathbf{C}) \end{array} \right) - \left( \begin{array}{c} \Re(Z_{T,SP})\\ \Im(Z_{T,SP}) \end{array} \right) \right) $$

The goal of the this section is to determine the elements of the Jacobian matrix $\mathbf{D}$.

Numerical Values of the Jacobian¶

Let's define a WEST ICRH antenna for a given frequency:

In [ ]:

f0 = 55  # MHz
freq0 = rf.Frequency(f0, f0, unit='MHz', npoints=1)
antenna = WestIcrhAntenna(frequency=freq0, 
                          front_face='../west_ic_antenna/data/Sparameters/front_faces/TOPICA/S_TSproto12_55MHz_Hmode_LAD6-2.5cm.s4p')
antenna.front_face_Rc()

Now we calculate the ideal matching point (for an half side only):

In [ ]:

C_match = antenna.match_one_side(f_match=f0*1e6, 
                                 side='left', solution_number=1)
print(C_match)

In [ ]:

power = [1,1e-6]
phase = [0,0]

Let's calculate the derivate of $Z_T$ around the match point:

In [ ]:

delta_C_tops = np.linspace(-10, 10, 51)
delta_C_bots = np.linspace(-10, 10, 51)

In [ ]:

# derivatives vs Ctop
z_Ts = []
z_T_SP = (2.9-0.3j)

for delta_C_top in delta_C_tops:
    # z_T
    _z_T = antenna.z_T(power, phase, Cs=[C_match[0]+delta_C_top, C_match[1], 150, 150])
    z_Ts.append(_z_T[:,0]) # left side
z_Ts = np.array(z_Ts).squeeze()
dz_Ts_dCtop = np.diff(z_Ts)

# normalized error and derivative
eps_Ctop = (z_T_SP - z_Ts)/z_Ts
deps_dCtop = np.diff(eps_Ctop)

# diff vs Cbot
z_Ts = []
for delta_C_bot in delta_C_bots:
    # z_T
    _z_T = antenna.z_T(power, phase, Cs=[C_match[0], C_match[1]+delta_C_bot, 150, 150])
    z_Ts.append(_z_T[:,0]) # left side
z_Ts = np.array(z_Ts).squeeze()/z_T_SP
dz_Ts_dCbot = np.diff(z_Ts)

# normalized error and derivative
eps_Cbot = (z_T_SP - z_Ts)/z_Ts
deps_dCbot = np.diff(eps_Cbot)

Below are the traces of the derivatives of $Z_T$ around the Set Point:

In [ ]:

fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)
axes[0,0].plot(delta_C_tops[:-1], dz_Ts_dCtop.real, label='$\partial \Re(z_T)/\partial C_{top}$')
axes[0,1].plot(delta_C_bots[:-1], dz_Ts_dCbot.real, label='$\partial \Re(z_T)/\partial C_{bot}$')
axes[1,0].plot(delta_C_tops[:-1], dz_Ts_dCtop.imag, label='$\partial \Im(z_T)/\partial C_{top}$')
axes[1,1].plot(delta_C_bots[:-1], dz_Ts_dCbot.imag, label='$\partial \Im(z_T)/\partial C_{bot}$')
axes[1,0].set_xlabel('$\delta C$')
axes[1,1].set_xlabel('$\delta C$')
[ax.legend() for ax in axes.ravel()]
[ax.axvline(0, color='gray', ls='--', alpha=0.8) for ax in axes.ravel()]
fig.subplots_adjust(wspace=0.0, hspace=0.0)

Below are the plot of the derivatives of the error signal $\varepsilon$ around the Set Point:

In [ ]:

fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)
axes[0,0].plot(delta_C_tops[:-1], deps_dCtop.real, label=r'$\partial \Re(\varepsilon)/\partial C_{top}$')
axes[0,1].plot(delta_C_bots[:-1], deps_dCbot.real, label=r'$\partial \Re(\varepsilon)/\partial C_{bot}$')
axes[1,0].plot(delta_C_tops[:-1], deps_dCtop.imag, label=r'$\partial \Im(\varepsilon)/\partial C_{top}$')
axes[1,1].plot(delta_C_bots[:-1], deps_dCbot.imag, label=r'$\partial \Im(\varepsilon)/\partial C_{bot}$')
axes[1,0].set_xlabel('$\delta C$')
axes[1,1].set_xlabel('$\delta C$')
[ax.legend() for ax in axes.ravel()]
[ax.axvline(0, color='gray', ls='--', alpha=0.8) for ax in axes.ravel()]
fig.subplots_adjust(wspace=0.0, hspace=0.0)

Below are the values of the derivative at the Set Point:

In [ ]:

idx = np.argmin(np.abs(delta_C_bots))
print('dz_T/dCtop:', dz_Ts_dCtop[idx])
print('dz_T/dCbot:', dz_Ts_dCbot[idx])

print('de/dCtop:', deps_dCtop[idx])
print('de/dCbot:', deps_dCbot[idx])

As the $\mathbf{D}$ matrix depends of the conditions (coupling, excitation), it is not possible to determine it exactly in advance in real experiments.

An approximative solution is to take only the sign of the derivatives: hence few iterations are needed until converging to the solution.

Note that this algorithm will diverge if the start point is not close enough from the optimal solution:

In [ ]:

D = np.array([[dz_Ts_dCtop[idx].real, dz_Ts_dCbot[idx].real],
              [dz_Ts_dCtop[idx].imag, dz_Ts_dCbot[idx].imag]])
print(np.sign(D))

So, for solution 1, one would obtain: $$ \left( \begin{array}{cc} -1 & -1 \\ -1 & +1 \end{array} \right) $$ and for solution 2: $$ \left( \begin{array}{cc} -1 & -1 \\ +1 & -1 \end{array} \right) $$

Note the sign change if one uses the Jacobian relative to the error signal $\mathbf{\varepsilon}$ instead of the Jacobian relative to the $z_T$.

Numerical Examples¶

In the example below, we setup and antenna and we look forward a solution by iterating on the capacitor predictor detailled below.

In [ ]:

f0 = 55  # MHz
freq0 = rf.Frequency(f0, f0, unit='MHz', npoints=1)
antenna = WestIcrhAntenna(frequency=freq0, 
                          front_face='../west_ic_antenna/data/Sparameters/front_faces/TOPICA/S_TSproto12_55MHz_Hmode_LAD6-2.5cm.s4p')
C_match_sol1 = antenna.match_one_side(f_match=f0*1e6, side='left', solution_number=1)
C_match_sol2 = antenna.match_one_side(f_match=f0*1e6, side='left', solution_number=2)
print(C_match_sol1)
print(C_match_sol2)

Matching a single side¶

In order to illustrate the convergence, we first calculate a capacitor map $\mathrm{VSWR}(C_{top}, C_{bot})$ for the left side around the ideal match point:

In [ ]:

delta_Cs = np.linspace(-10, +10, num=31, endpoint=True)

In [ ]:

C_top_lefts, C_bot_lefts = np.meshgrid(C_match_sol1[0] + delta_Cs, 
                                       C_match_sol1[1] + delta_Cs)
print(C_top_lefts.size)

In [ ]:

power = [1, 1e-12]  # small value on right side to avoid division by zero
phase = [0, 0]

vswrs = []
for (C_top, C_bot) in tqdm(np.nditer([C_top_lefts, C_bot_lefts]), total=C_top_lefts.size):
    _vswr = antenna.vswr_act(power, phase, Cs=[C_top, C_bot, 120, 120])
    vswrs.append(_vswr)
vswrs = np.array(vswrs).squeeze()
vswrs_left = vswrs[:,0].reshape(C_top_lefts.shape)

Which looks like that:

In [ ]:

fig, ax = plt.subplots()
cs=ax.contourf(C_bot_lefts, C_top_lefts, vswrs_left, np.linspace(1, 6, 41))
cs2=ax.contour(C_bot_lefts, C_top_lefts, vswrs_left, np.linspace(1, 6, 11), colors='k', alpha=0.6)
ax.clabel(cs2, inline=1, fontsize=10)
ax.set_xlabel('$C_{bot}$ [pF]')
ax.set_ylabel('$C_{top}$ [pF]')
ax.set_title('SWR at feeding line - left side')
ax.grid(True, alpha=0.2)
ax.plot(C_match_sol1[1], C_match_sol1[0], 'o', color='C1', label="solution 1 ($C_{top} > C_{bot}$)")
ax.plot(C_match_sol2[1], C_match_sol2[0], 'o', color='C2', label="solution 2 ($C_{top} > C_{bot}$)")
ax.legend()
fig.colorbar(cs)

Now let's define a starting point and iterate on the capacitor predictor. Playing with the code below, you'll see rapidly that the convergence is heavily linked to the starting point. Using a starting point far from the expected solution or with capacitor inversed with respected to the desired solution will lead to the divergence of the algorithm:

In [ ]:

C0_start = [60,40,120,120]
sol_num = 1

C0 = list(C0_start)  # new list to avoid reference
C0s = []
cont = True
iterations = 0 
while cont:
    C_next_left, C_next_right, eps = antenna.capacitor_predictor(power, phase, Cs=C0, solution_number=sol_num)
    C0 = [C_next_left.squeeze()[0], C_next_left.squeeze()[1], 120, 120]

    iterations += 1
    if iterations > 1 and (np.abs(C0s[-1][0] - C0[0]) < 0.01) and (np.abs(C0s[-1][1] - C0[1]) < 0.01):
        cont = False
    if iterations > 60:
        cont = False
    C0s.append([C0[0], C0[1]])
print(f'Stopped after {iterations} iterations')    

Let's illustrate the trace of the matching convergence:

In [ ]:

fig, ax = plt.subplots()
cs=ax.contourf(C_bot_lefts, C_top_lefts, vswrs_left, np.linspace(1, 6, 41))
cs2=ax.contour(C_bot_lefts, C_top_lefts, vswrs_left, np.linspace(1, 6, 11), colors='k', alpha=0.6)
ax.clabel(cs2, inline=1, fontsize=10)
ax.set_xlabel('$C_{bot}$ [pF]')
ax.set_ylabel('$C_{top}$ [pF]')
ax.set_title('SWR at feeding line - left side')
ax.grid(True, alpha=0.2)
ax.plot(C_match_sol1[1], C_match_sol1[0], 'o', color='C1', label="solution 1 ($C_{top} > C_{bot}$)")
ax.plot(C_match_sol2[1], C_match_sol2[0], 'o', color='C2', label="solution 2 ($C_{top} > C_{bot}$)")
ax.legend()
fig.colorbar(cs)


C0s = np.array(C0s)
ax.plot(C0_start[1], C0_start[0], 'x', color='r')
ax.plot(C0s[:,1], C0s[:,0], '-x', color='r', alpha=0.8)
print(C0s[-1])

Below is another representation of the convergence:

In [ ]:

fig, ax = plt.subplots()
ax.plot(np.arange(len(C0s)), C0s, lw=2)
ax.set_xlabel('# Steps')
ax.set_ylabel('Capacitance [pF]')
ax.grid(True, alpha=.4)
ax.set_title('Capacitor Matching Alg Convergence')
ax.axhline(eval(f'C_match_sol{sol_num}')[0], ls='--', color='gray')
ax.axhline(eval(f'C_match_sol{sol_num}')[1], ls='--', color='gray')

Matching both sides at the same time¶

In [ ]:

power = [1, 1]
phase = [0, np.pi]
# reference solution we wish to obtain
C_opt_dipole = antenna.match_both_sides(f_match=55e6, power=power, phase=phase)

In [ ]:

C0_start = [60, 40, 60, 40] # note that the start point matches solution 1 (Ctop>Cbot)
sol_num = 1

C0 = list(C0_start)  # new list to avoid reference
C0s = []
cont = True
iterations = 0 
while cont:
    C_next_left, C_next_right, eps = antenna.capacitor_predictor(power, phase, Cs=C0, solution_number=sol_num)
    C0 = [C_next_left.squeeze()[0], C_next_left.squeeze()[1], 
          C_next_right.squeeze()[0], C_next_right.squeeze()[1]]
    print(C0)
    iterations += 1
    if iterations > 1 and (np.abs(C0s[-1][0] - C0[0]) < 0.01) and (np.abs(C0s[-1][1] - C0[1]) < 0.01):
        cont = False
    if iterations > 60:
        cont = False
    C0s.append([C0[0], C0[1], C0[2], C0[3]])
print(f'Stopped after {iterations} iterations')    

In [ ]:

fig, ax = plt.subplots()
ax.plot(np.arange(len(C0s)), C0s, lw=2)
ax.set_xlabel('# Steps')
ax.set_ylabel('Capacitance [pF]')
ax.grid(True, alpha=.4)
ax.set_title('Capacitor Matching Alg Convergence')
[ax.axhline(C, ls='--', color=f'C{idx}') for idx, C in enumerate(C_opt_dipole)]

The antenna is matched:

In [ ]:

antenna.vswr_act(power, phase, Cs=C0s[-1])

The automatic matching algorithm is implemented in the matching_both_sides_iterative method of the antenna object:

In [ ]:

C_opt_2 = antenna.match_both_sides_iterative(f_match=55e6, power=power, phase=phase, Cs=[50, 50, 50, 50])

In [ ]:

antenna.vswr_act(power, phase, Cs=C_opt_2)

Correspondance Between Error Signal Amplitude and Capacitance¶

In the control room, the IC operator job is often to improve the antenna tuning for the next shot. Because turning on the matching feedback is not always desirable, another option is to look to the error signals $\mathbf{\varepsilon}$ generated in the latest pulse and to deduce from them the correction to apply to the capacitors.

First, let's build a half-matched antenna:

In [ ]:

f0 = 55  # MHz
freq0 = rf.Frequency(f0, f0, unit='MHz', npoints=1)
antenna = WestIcrhAntenna(frequency=freq0, 
                          front_face='../west_ic_antenna/data/Sparameters/front_faces/TOPICA/S_TSproto12_55MHz_Hmode_LAD6-2.5cm.s4p')
C_match = antenna.match_one_side(f_match=f0*1e6, 
                                 side='left', solution_number=1)
print(C_match)

Now, let's depart one capacitor from the ideal point and look to the evolution of the error signal:

In [ ]:

delta_Cs = np.linspace(-10, +10, 11)

_Cs = []
for delta_C in tqdm(delta_Cs):
    Cs = [C_match[0]+delta_C, *C_match[1:]]
    _C = antenna.capacitor_predictor(power, phase, Cs)[0].squeeze()
    _Cs.append(_C)
_Cs = np.array(_Cs)

In [ ]:

fig, ax = plt.subplots()
ax.plot(delta_Cs, _Cs-C_match[:2])
ax.plot([-10, 10], [-10, 10])

Velocity Controler¶

The velocity controler matrix $\mathbf{T}$ is in fact $\mathbf{D}^{-1}$, since the velocity $\mathbf{V}$ is proportional to the change $\delta C$:

$$ \delta C = \mathbf{D}^{-1} (z_T - z_{T,SP}) $$

so that $$ \mathbf{V} = k \times \delta C = k \mathbf{D}^{-1} (z_T - z_{T,SP}) = \mathbf{T} (z_T - z_{T,SP}) $$ Following the previous result on $\mathbf{D}$, the matrix $\mathbf{T}$ can be parametrized as:

$$ \mathbf{T} = k\times \left( \begin{array}{cc} -S \cos \alpha & S K_{ri} \sin\alpha \\ -\sin\alpha & - K_{ri} \cos\alpha \end{array} \right) $$

Where $k$ is gain, $S$ and $\alpha$ two parameters used to define the choice of the matching solution ("1" for $C_{top}>C_{bot}$, "2" for the opposite). $K_{ri}$ is used to define the relative weight between real and imaginary part of the error signal $\mathbf{\varepsilon}$.

In [ ]:

from IPython.core.display import HTML
def _set_css_style(css_file_path):
    """
    Read the custom CSS file and load it into Jupyter
    Pass the file path to the CSS file
    """
    styles = open(css_file_path, "r").read()
    s = '<style>%s</style>' % styles
    return HTML(s)

_set_css_style('custom.css')

In [ ]: