Model-Based Optimization of the Field-Null Configuration for Robust Plasma Breakdown on the HL-3 Tokamak

Abstract

This paper introduces a self-consistent field-null optimization algorithm of a poloidal magnetic field that precisely accounts for the influence of vacuum vessel eddy currents. Building on existing poloidal field (PF) coil currents, the algorithm can refine these waveforms to achieve various target field-null configurations. Firstly, based on the TokSys toolbox, a response model, including the PF coils and vacuum vessel circuits for the HL-3 tokamak, is developed under the MATLAB^® and Simulink^™ framework. The resistivity parameters of the model are calibrated using experimental data obtained from single-coil discharge tests. Subsequently, an iterative method was employed to simultaneously solve the dynamic field-null optimization problem within a specified spatial region and precisely account for the effect of passive eddy currents. Typically, $\left|B_{\perp }\right|\le 1 G$ within a large area can be obtained with this iterative scheme, which can be stably sustained for over 15 milliseconds to ensure the robustness of breakdown. Finally, a low-pass filtered PID controller is applied to the model to achieve precise control of the PF coils currents, confirming the feasibility of implementing the proposed algorithm in real experiments.

1. Introduction

Nuclear fusion can generate tremendous energy with a very small amount of fuel and is often considered the ultimate energy source for humanity. Up to now, tokamak is considered to be the most promising way to produce fusion energy, which is a doughnut shaped device with a strong magnetic field inside to confine the fusion plasmas with an extremely high temperature [1]. In tokamaks, fusion plasmas are usually initiated by breaking down the neutral gases through a Townsend avalanching process by applying an inductive toroidal electric field in a region with a low poloidal magnetic field [1,2]. It is essential in tokamaks to optimize the breakdown process, not only because it is sensitive to various factors and requires careful design but also due to its close link with the later stages of plasma initialization, such as impurity burn-through, plasma equilibrium formation, current ramp-up and shape control [3]. Most of the major tokamaks around the world have made extensive efforts in improving the reliability and quality of plasma breakdown, like JET, DIII-D, EAST and ITER, etc. [4,5,6,7]. Generally, the field-null configuration is a critical issue in breakdown optimization, as the breakdown condition of tokamaks is written as [8]:

$$E_{min}\left(\mathrm{V}/\mathrm{m}\right)={\displaystyle \frac{1.25\times {10}^4 p \left(Torr\right)}{\ln \left[510 p\left(Torr\right)L\left(\mathrm{m}\right)\right]}}$$

(1)

where $L={0.25a}_{eff}{B}_T/B_{\perp }$ is the effective connection length. This equation shows that obtaining a field-null configuration, i.e., minimizing the poloidal magnetic field $B_{\perp }$, can significantly reduce the required toroidal electric field for a successful plasma breakdown, which is essential as the available toroidal electric field is quite limited on large super-conducting machines, such as ITER [9]. Optimizing the field-null configuration is delicate work, as $B_{\perp }$ is determined not only by currents in the active controlled magnetic coils but also by the eddy currents in various passive conductive structures. Special attention should be paid to the treatment of passive eddy currents, which could be as large as tens to hundreds of kilo-amperes and can significantly influence $B_{\perp }$ during the breakdown stage due to the high inductive voltage, as has been demonstrated in the WEST tokamak [10].

To address this problem, various methods have been proposed. In the early 1990 s, a model-based approach was developed to optimize the coil voltage waveforms for tokamak startup [11,12]. It utilizes the eigenvalue transformation [13] to express the field values like $B_{\perp }$ linearly in terms of coil voltages, with the effect of eddy currents implicitly included. This transfers the optimization problem of $B_{\perp }$ to minimizing a quadratic cost function, which is much easier to solve. This approach has been widely used for designing field-null configurations on several devices, like ITER [12], EAST, KSTAR [14] and HT-6M [15]. Based on similar principles, more integrated approaches have also been proposed to optimize the whole startup scenario including breakdown and the following processes, like plasma current ramp-up and equilibrium formation [7,9,16,17]. Yet, these approaches rely on the assumption of piecewise linear coil voltage waveforms, which could pose a limit for both the tuning flexibility and the consistence with arbitrary waveforms when applying feedback-control algorithms. Other approaches include the one based on flux expansion in DIII-D [5] and adjusting the currents of a specific set of coils empirically to compensate for the influence of eddy currents, such as on HL-3 (previously named HL-2M) [18,19]. Neither of these approaches has an explicit, self-consistent treatment of eddy currents, which limits the accuracy of $B_{\perp }$ compensation and the quality of designed field-null configurations. On the other hand, issues related to feedback control often lack mentioning in the literature, which requires attention as it is fundamental in filling the gap between modeling and experiments, and the deviation from feedforward trajectories could significantly impact the field-null configuration.

In this paper, a model-based scheme for optimizing the field-null configuration is proposed and applied on HL-3, which is a medium-sized copper tokamak aimed at exploring the route to high-performance plasma operation [20,21]. The scheme is based on an integrated forward simulation model of the tokamak circuit system, including both the active coils and passive conductors, which is developed within the MATLAB^® Simulink™ framework. The scheme does not require additional assumptions of coil voltage waveforms, such as being piecewise linear, which could make it more straightforward and flexible when designing the coil current and voltage waveforms. Similar to other model-based schemes, this one has the advantage of being self-consistent as well, and the passive eddy currents are precisely treated by applying an iterative optimization approach. A large field-null area with $\left|B_{\perp }\right|\le 1 G$ is obtained when applying this scheme for a breakdown scenario design on HL-3, which can be steadily sustained for over 15 milliseconds to guarantee the robustness of breakdown. The corresponding coil current waveforms are also provided as the control targets to be used in experiments. To verify the feasibility of the proposed scheme in closed-loop environments, a feedback controller based on a low-pass filtered PID algorithm is integrated with the forward model to perform simulation tests, and the result confirms the high potential of applying this scheme in HL-3 experiments. In addition, for newly-built devices, like ITER, in the future, the proposed scheme could also provide an alternative and flexible approach in breakdown scenario design for obtaining robust tokamak startups.

The remainder of this paper is organized as follows. Section 2 introduces the forward model of coils and passive conductive structures on HL-3 based on the RL circuit principle and its validation results. Section 3 proposes an iterative scheme for field-null optimization based on the model in Section 2 and presents some preliminary designs of field-null configurations for various purposes on HL-3. The corresponding coil current waveforms are also provided as the target for feedback control in not only simulations but also actual experiments. Section 4 deals with the feedback controller and builds a closed simulation control loop based on a filtered PID algorithm, which is then used to perform closed-loop control simulations that verify the feasibility of designed coil current waveforms. Finally, Section 5 presents the conclusion and some discussions.

2. Modeling of the Poloidal Field Coil System on HL-3

The poloidal field system of the HL-3 tokamak consists of 16 PF coils and a central solenoid (CS) coil [22]. The cross-sectional diagram of the coils and vacuum vessel is shown in Figure 1. The 16 PF coils are located inside the toroidal field coils (TF coils) and outside the vacuum vessel, symmetrically distributed along the midplane of the HL-3. To facilitate control, the PF and CS coils are powered by separate power supplies [23].

2.1. Mathematical Model

2.1.1. Conductors Model

The design and optimization of the field-null configuration primarily focus on the poloidal component of the magnetic field. The poloidal magnetic field component induced by the TF coils is negligible according to a careful analysis of the magnetic probe data. Consequently, only the CS and PF coils and passive conductor structures are needed for the modeling. Primarily, the passive conductor structures include only the vacuum vessel, as the first wall modules on HL-3 are isolated from each other, and no passive coils are present inside the vessel.

The HL-3 tokamak is toroidally axisymmetric, where all CS coils and PF coils can be approximately treated as multilayer, tightly wound circular coils with rectangular or parallelogram cross-sections, distributed in parallel and coaxially. Assuming that the vacuum vessel material has a uniform medium, it can be equivalently modeled as a number of closed circular passive conductor loops distributed coaxially in the poloidal direction and connected in parallel. In the following model, the number of passive conductor loops is chosen to be 88 to obtain high accuracy while retaining a reasonable computational cost.

After simplifying the conductors circuit, the relationship between the voltages and currents in the conductors can be expressed as follows:

$$V_i=R_i{I}_i+L_i{\displaystyle \frac{d{I}_i}{dt}}+{\sum }_{j=1}^{j=105}{M}_{ij}{\displaystyle \frac{d{I}_j}{dt}} \left(i \ne j\right)$$

(2)

The subscript $i=1, 2,\dots ,105$ represents all the CS, PF coils and passive conductor loops, among which $i=1$ represents the CS coil, $i=2~17$ represents the 16 PF coils, and $i=18~105$ represents the 88 passive conductors. Here, $I_i$ and $V_i$ are the current and terminal voltage of the $i$th conductor, respectively, $R_i$ is the total resistance of the $i$th coil and its transmission line, $L_i$ is the self-inductance of the $i$th coil, and $M_{ij}$ is the mutual inductance between the $i$th and $j$th coil. Apparently, $V_i$ is zero for $i=18~105$ as they are passive conductors and there is no power supply that energizes them.

Equation (1) can be rewritten in its matrix form:

$$\left[\begin{array}{c}\begin{array}{c}{V}_1\\ V_2\\ ⋮\end{array}\\ V_{105}\end{array}\right]=\left[\begin{array}{cccc}{R}_1& 0& \dots & 0\\ 0& R_2& \dots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \dots & R_{105}\end{array}\right]\left[\begin{array}{c}{I}_1\\ I_2\\ ⋮\\ I_{105}\end{array}\right]+\left[\begin{array}{cccc}{M}_{1,1}& M_{1,2}& \dots & M_{1,105}\\ M_{2,1}& M_{2,2}& \dots & M_{2,105}\\ ⋮& ⋮& \ddots & ⋮\\ M_{105,1}& M_{105,2}& \dots & M_{105,105}\end{array}\right]\left[\begin{array}{c}{\dot{I}}_1\\ {\dot{I}}_2\\ ⋮\\ {\dot{I}}_{105}\end{array}\right]$$

(3)

which can be simplified as follows:

$$\mathit{V}=\mathbf{R}\mathit{I}+\mathbf{M}\dot{\mathit{I}}$$

(4)

In Equation (4), $\mathit{I}$ and $\mathit{V}$ are the$105\times 1$current and voltage vector, respectively, and $\dot{\mathit{I}}\equiv d\mathit{I}/dt$ is the time derivative of coil currents $\mathit{I}$, $\mathbf{M}$ is a $105\times 105$ symmetric matrix containing all of the self-and-mutual inductance values, and $\mathbf{R}$ is a $105\times 105$diagonal matrix containing all of the resistance values of coils and passive conductor loops.

The Equation (4) can be rewritten as follows:

$$\dot{\mathit{I}}=\mathbf{A}\mathit{I}+\mathbf{B}\mathit{V}$$

(5)

In Equation (5), $\mathbf{A}=-{\mathbf{M}}^{-1}\mathbf{R}, \mathbf{B}={\mathbf{M}}^{-1}$.

The tool used to calculate the self-inductance and mutual inductance of the coils is TokSys, which calculates the self-inductance of each coil and the mutual inductance between coils based on the Biot-Savart law and elliptic integrals [24].

2.1.2. Power Supply Model

The CS/PF power supply system utilizes a four-quadrant controlled thyristor rectifier bridge structure with an internal circulating current to energize the coils. The triggering angle of the thyristors is controlled by the Plasma Control System (PCS) to adjust the terminal voltage of the PF coils. In the control model, the power supply can be regarded as an inertial element with delay [25]. The input is the output voltage of the power supply requested by the PCS, and the output is the actual output voltage value. Assuming that the proportional coefficient is 1, the response time constant is ${\tau }_{ps}$, and the power supply delay time is $\tau$, and the transfer function of the power supply system can be expressed as follows [25]:

$$G\left(s\right)=\frac{e^{-\tau s}}{{s\tau }_{ps}+1}$$

(6)

2.2. Tuning of the Model Parameters

Considering the fact that PF coils are energized at different times in a discharge, Equation (5) is treated as a time-varying state equation with parameter changes. In addition, the accuracy of the model is the foundation for performing zero-field optimization and developing subsequent control algorithms. The parameters of the model calculated by TokSys do not take into account the errors in the actual construction process, and there are some differences with the actual situation. To make sure that the simulation model fits the hardware characteristics of the actual system, the mutual inductance and resistance coefficients are calibrated using single-coil discharge experimental data. The input is the actual voltage supplied during the experiment, while the outputs are the coil current and the eddy current in the vacuum chamber. The voltage/current of the flat top section is used to calculate the resistance of the coil, and then, the self-inductance is adjusted to ensure that the simulation output is consistent with the experimental results. The calibrated model was then validated against a normal plasma discharge shot, and the results are shown in Figure 2. It can be seen that both in the current ramp-up phase and the flat-top phase, the corrected response model’s output closely fits the experimental data with an error level less than 5%. This demonstrates the reliability of the time-varying parameter model, which forms a solid basis for the following field-null optimization work.

3. Optimization of the Field-Null Configuration

Robust plasma breakdown in tokamaks relies significantly on the quality of the field-null configuration, i.e., the absolute value of the poloidal field $\left|B_{\perp }\right|$, the area and shape of the field-null region, and its alignment with the high-loop-voltage region. In particular, robust breakdown with a low toroidal electric field $E_{\varphi }\le 0.3 \mathrm{V}/\mathrm{m}$ is desirable on large tokamaks like ITER [9], due to the limited change rate of coil currents for super conductors. However, for now, on HL-3, the breakdown is not very reliable even with $\left|E_{\varphi }\right|$ as high as $1.0~1.2 \mathrm{V}/\mathrm{m}$. Thus, there is an urgent need for improving the breakdown scenario on HL-3. One of the key issues in breakdown improvement is how to compensate for the influence of eddy currents in the passive conductors to obtain a good field-null configuration. Based on the simulation model in Section 2, an iterative scheme is introduced to minimize the average $\left|B_{\perp }\right|$ in a given region and simultaneously provide the corresponding PF coil current waveforms. The scheme is simple enough to be calculated swiftly, and the eddy currents in passive structures are self-consistently treated to ensure a reliable solution for compensation. The remaining part of this section tells the details of the optimization scheme and typical examples of optimized field-null configurations.

3.1. The Iterative Optimizing Scheme

The optimization problem of the field-null configuration could be abstracted to minimize the averaged $\left|B_{\perp }^2\right|$ inside a given region. Expressed more explicitly, it is to minimize an optimizing function $f$ by discretizing the tokamak cross-section into a set of grids as

$$f\left(\mathit{I}\right)={\sum }_j^{j\in \mathbf{A}}\left|B_{p,j}^2\left(\mathit{I}\right)\right|={\sum }_j^{j\in \mathbf{A}}|B_{R,j}^2\left(\mathit{I}\right)+B_{Z,j}^2\left(\mathit{I}\right)|$$

(7)

Here, $\mathit{I}$ is the currents in coils and passive conductors, which is the same as in Equation (5), $\mathbf{A}$ represents the range of the desired field-null region, the subscript $j$ refers to the grid number and $B_{R,j}, B_{z,j}$and $B_{p,j}$are the horizontal component, vertical component and the sum of the poloidal magnetic field on the $j$th grid, respectively. By utilizing Green’s function, Equation (7) can be rewritten as a multi-variable, quadratic function:

$$f\left(\mathit{I}\right)={\sum }_j^{j\in \mathbf{A}}|B_{R,j}^2\left(\mathit{I}\right)+B_{Z,j}^2\left(\mathit{I}\right)|={\sum }_j^{j\in \mathbf{A}}\left\{{\left[{\sum }_{i=1}^{105}{G}_{R,ij}{I}_i\right]}^2+{\left[{\sum }_{i=1}^{105}{G}_{Z,ij}{I}_i\right]}^2\right\}$$

(8)

Here, $G_{R,ij}$and$G_{Z,ij}$are Green’s functions that represent the $R$ and $Z$ components of the poloidal field induced by a 1-ampere current in the $i$th conductor on the $j$th grid, and $I_i$ is the current in the $i$th conductor.

Various methods could be applied to solve the minimization problem of $f\left(\mathit{I}\right)$ as expressed by Equation (8). For simplicity, taking the quadratic form of $f$ into consideration, a direct scheme is to calculate the solution of ${\displaystyle \frac{\partial f}{\partial \mathit{I}}}=0$ that simultaneously satisfies the positive-definite condition of the Hessian matrix, e.g., $\left|{\displaystyle \frac{{\partial }^{2}f}{\partial I_{i1}\partial I_{i2}}}\right|>0$. The problem ${\displaystyle \frac{\partial f}{\partial \mathit{I}}}=0$ could be written as a set of linear equations:

$$\mathbf{P}\cdot \mathit{I}=0$$

(9)

Here, $\mathbf{P}\equiv P_{i^{\prime }i}$ is an $105\times 105$ matrix with $P_{i^{\prime }i}=2{\sum }_j^{j\in \mathbf{A}}\left(G_{R, i^{\prime }j}\cdot G_{R,ij}+G_{Z,i^{\prime }j}\cdot G_{Z,ij}\right)$. Equation (9) is homogeneous and has only a trivial solution $I_i=0$ if all of the $I_i$ values are unknowns. However, an exclusive, non-trivial solution usually exists if some elements in $\mathit{I}$ are treated as knowns. Without the loss of generality, assume that $\mathit{I}$ could be rewritten into the combination of two sub-vectors as

$$\mathit{I}=\left[\begin{array}{c}{\mathit{I}}_{\mathit{u}\mathit{n}\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}\\ {\mathit{I}}_{\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}\end{array}\right]$$

(10)

where ${\mathit{I}}_{\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}$ and ${\mathit{I}}_{\mathit{u}\mathit{n}\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}$ refer to the sub-vectors containing the known and unknown elements of vector $\mathit{I}$, respectively. Then, Equation (9) can be rewritten as

$${\mathbf{P}}_{\mathbf{1}}\cdot {\mathit{I}}_{\mathit{u}\mathit{n}\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}={\mathbf{P}}_{\mathbf{2}}\cdot {\mathit{I}}_{\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}$$

(11)

where ${\mathbf{P}}_{\mathbf{1}}\equiv P_{unknown,unknown}$ is a square, symmetric sub-matrix of $\mathbf{P}$ that extracts the elements in $\mathbf{P}$ with both the same row and column subscripts as ${\mathit{I}}_{\mathit{u}\mathit{n}\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}$ within $\mathit{I}$, and ${\mathbf{P}}_2\equiv P_{unknown,known}$ is another sub-matrix of $\mathbf{P}$ that extracts the elements in $\mathbf{P}$ with the same row and columns subscripts as ${\mathit{I}}_{\mathit{u}\mathit{n}\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}$ and ${\mathit{I}}_{\mathit{k}\mathit{n}\mathit{o}\mathit{w}\mathit{n}}$ within $\mathit{I}$, respectively.

Equation (11) is useful when some conductor currents can be treated as knowns, i.e., when the CS current waveform is fixed to optimize the field-null configuration with a certain level of $E_{\varphi }$, and setting some PF currents as zero to reduce the complexity of the optimization problem. Under this circumstance, the CS current $I_{CS}$ and part of the PF currents $I_{PF,known}$ can be treated as known, while the remaining part of PF currents $I_{PF,unknown}$ is to be solved. The eddy currents in passive conductors could pose a problem in solving this problem, as these currents are determined by Equation (5) in the form of time derivatives that interfere with the linear nature of Equation (11). Here, a simple scheme is introduced to also treat the eddy currents as “known” and to obtain them by iteratively solving Equations (11) and (5), as shown by the flow chart in Figure 3:

The scheme is proven to be efficient and usually converges in roughly 30 iterations. Figure 4 shows an example of minimizing the $\left|B_{\perp }\right|$ in a circular region with $R=1.78 \mathrm{m}$ and $r=0.65 \mathrm{m}$, where Figure 4a–c shows the contour plot of optimized $\left|B_{\perp }\right|$ in iteration cycles 1, 29 and 30, respectively. Figure 4d shows the sum of eddy currents in iteration cycles 1, 29 and 30. It can be seen that the difference between $I_{Eddy}$ in iteration cycle 1 and 29 is large, but the difference between cycles 29 and 30 is negligible, and the $\left|B_{\perp }\right|$ distribution is almost the same between cycles 29 and 30, indicating that this scheme converges fairly well. Thus, the calculation of eddy currents could be regarded as self-consistent.

3.2. Typical Examples of Field-Null Optimization

Section 3.1 presents a scheme to obtain the PF coil current waveforms in order to obtain an optimum field-null configuration, which is valid for both the circumstances of static and evolving coil currents with a self-consistent eddy current in passive conductors. This provides a useful tool for designing not only static but also evolving field-null configurations, which is particularly useful when the eddy current is high and requires careful compensation. In addition, by choosing different regions $\mathbf{A}$ in Equations (7)–(11), one can conveniently design field-null configurations for various purposes.

For example, if the field-null area is desired to be as large as possible to ensure robust breakdown, one can choose $\mathbf{A}$ to be a big circular region as shown by the red dashed circle line in Figure 5a. The CS current is fixed to be the control target waveform used in HL-3 experiments, which is treated as known, and PF6 + PF8 coils are used to form the field-null configuration for simplicity. By following the workflow shown in Figure 3, a set of PF6 and PF8 current waveforms is obtained for minimizing the averaged $\left|B_{\perp }\right|$ inside region $\mathbf{A}$, as shown in Figure 5b. These waveforms, together with the preset CS current waveform, can be used as control targets in actual experiments. In another case, the PF4 and PF7 coils are used additionally as unknowns together with PF6 and PF8, a much better field-null configuration is obtained, as shown in Figure 5c. The corresponding CS and PF current waveforms are presented in Figure 5d.

In other situations, there might be requirements to prevent the plasma from losing vertical stability at the early stage of current ramp-up, and the oblate elliptic region as shown in Figure 6a could be chosen to make the plasma breakdown primarily near the midplane, avoiding the formation of an early plasma shape that is too elongated and could be vertically unstable. For some harsh circumstances, like the wall condition being bad and the impurity level being high, a field-null region shown in Figure 6c could be selected to align with the high-loop-voltage region on the high field side (HFS) of a tokamak. The CS and PF currents corresponding to Figure 6a,c are shown in Figure 6b,d, respectively. These waveforms can provide the basis for feedback-control simulations and could greatly reduce the time cost of tuning the field-null configuration by modifying coil currents only experimentally.

4. Feedback Control System and Field-Null Implementation

Section 3 introduces the scheme for field-null optimization and gives the control target waveforms of CS and PF currents for various field-null configurations. To realize these currents in real experiments, precise feedback (FB) control must be used. In this section, a FB controller based on the Low-Pass Filtered PID principle [26,27] is introduced and used for simulating the controlling process of CS and PF coil currents. The result shows that the designed target coil current waveforms can be accurately tracked by the controller, which further validates the feasibility of these waveforms as control targets to be used in actual HL-3 experiments.

4.1. Feedback Control System

PID [28] is a widely used feedback-control method in industrial systems. On HL-3, PID controllers are integrated into the Plasma Control System (PCS) to handle different control objects. In actual plasma discharges, noises in diagnostic signals are inevitable and usually can lower the performance of PID controllers. To solve this problem, a low-pass filter can be added before the signal goes through PID controllers. In this way, the noises in the diagnostic signals can be eliminated, which improves the performance of the controller, and then, the processed error signals are proportionally multiplied, differentiated and integrated to form the P, D and I control quantities, respectively. The block diagram of the feedback control algorithm is shown in Figure 7.

Considering that the control period of PCS on HL-3 is 1 millisecond, and discrete signals are used in actual control, the filtering algorithms employed in this control algorithm are all digital filtering algorithms. During the calculation process, the impulse transfer function $G\left(s\right)$ of the analog filter can first be designed based on technical specifications, and then, the system function $H\left(z\right)$ of the corresponding discrete system can be obtained by applying the bilinear transform method.

4.1.1. Low-Pass Filtered PID Principle

Assuming that the time constant of an analog low-pass filter is ${\tau }_1=R_1{C}_1$, its impulse transfer function $G_1\left(s\right)$ is given by the following:

$$G_1\left(s\right)={\displaystyle \frac{1}{1+{\tau }_{1}s}}$$

(12)

The mapping from the $s$-plane to the $z$-plane is given by $s={\displaystyle \frac{2}{∆t}}{\displaystyle \frac{z-1}{z+1}}$, where $∆t$ is the sampling interval. The discrete system function $H_1\left(z\right)$ can be derived as follows:

$$H_1\left(z\right)={\displaystyle \frac{Y\left(z\right)}{U\left(z\right)}}={\displaystyle \frac{L_1(1+z^{-1})}{1+{L_{2}z}^{-1}}}$$

(13)

Taking the inverse $z$-transform of Equation (13), the discrete-time difference equation can be derived as follows:

$$Y\left(t\right)=L_{1}U\left(t\right)+L_{1}U\left(t-\mathsf{\Delta }t\right)-L_{2}Y\left(t-\mathsf{\Delta }t\right)$$

(14)

where $L_1={\displaystyle \frac{{\tau }_1^{*}}{{\tau }_1^{*}+1}}$, $L_2={\displaystyle \frac{{\tau }_1^{*}-1}{{\tau }_1^{*}+1}}$,${\tau }_1^{*}={\displaystyle \frac{\mathsf{\Delta }t}{2{\tau }_1}}$, $U\left(t\right)$ and $Y\left(t\right)$ represent the input and output signals at time $t$. Note that Equation (14) can be applied to both the low-pass filters of vectors $\mathbf{e}$ and $\mathbf{v}$ in Figure 7.

Then let us consider the high-pass filter. By applying the same method as that used for deriving Equation (14), one can obtain the discrete-time difference equation of a high-pass filter as

$$Y\left(t\right)=D_{1}U\left(t\right)-D_{1}U\left(t-∆t\right)-D_{2}Y\left(t-∆t\right)$$

(15)

where ${\tau }_2=R_2{C}_2$ is the time constant of an analog high-pass filter, $D_1={\displaystyle \frac{1}{{\tau }_2^{*}+1}}$, $D_2={\displaystyle \frac{{\tau }_2^{*}-1}{{\tau }_2^{*}+1}}$, ${\tau }_2^{*}={\displaystyle \frac{\mathsf{\Delta }t}{2{\tau }_2}}$, and $U$, Y has the same meaning as in Equation (14).

Using Equation (14) and (15), the I/O relationship of the low-pass filtered PID control algorithm can be written as follows:

$$\mathbf{v}\left(t\right)=L_1\mathbf{e}\left(t\right)+L_1\mathbf{e}\left(t-\Delta t\right)-L_2\mathbf{v}\left(t-\Delta t\right)$$

$$\mathbf{d}\left(t\right)=D_1\mathbf{v}\left(t\right)-D_1\mathbf{v}\left(t-∆t\right)-D_2\mathbf{d}\left(t-∆t\right)$$

$$\mathbf{i}\left(t\right)=I_1\mathbf{v}\left(t\right)+I_1\mathbf{v}\left(t-∆t\right)-I_2 \mathbf{i}\left(t-∆t\right)$$

$$u=K_P\ast \mathbf{v}\left(t\right)+K_D\ast \mathbf{d}\left(t\right)+K_I\ast \mathbf{i}\left(t\right)$$

(16)

Here, the expressions for $L_1$, $L_2$, $I_1$ and $I_2$ follow the form of the low-pass filter coefficients in Equation (14), $D_1$, $D_2$ follows the ones in Equation (15) and $K_p$, $K_D$ and $K_I$ are the gains of P, D and I, respectively.

4.1.2. Closed-Loop Control System in Simulink

The closed-loop control system for the CS and PF coils is built using Simulink^™ (version: 24.1), with the state-space equation set (2), which were tuned in Section 2.2, serving as the model for the coils. The Varying State Space module in Simulink is used to process the time-dependent parameters of the state-space equations into a structure. The From Workspace module is used to import data from the MATLAB (version: 24.1) workspace. The closed-loop control system structure, taking the CS coil as an example, is shown in Figure 8.

4.2. Simulation Results

By utilizing the closed-loop control system, simulation experiments on coil current FB control can be performed to determine whether the desired field-null configuration can be readily obtained, which form the last step before putting the designed field-null recipes into real HL-3 experiments. A typical example of these simulation experiments is shown in Figure 9, in which the field-null recipe shown by Figure 6a,b is selected as the control target, as such a configuration has not been realized in real experiments. The dashed lines in Figure 9a represent the simulatively controlled coil currents from the feedback control model, while the solid lines represent the target currents obtained by solving the field-null optimization problem, as shown in Figure 6b. It can be seen that the output values of coil currents from the low-pass-filtered PID controller can perfectly match the target current. As shown by Figure 9b, the typical control error is mostly below 0.3 kA for CS, 0.05 kA for PF4 and PF6 and 0.005 kA for PF7 and PF8 coils. This accuracy ensures that the subsequent $\left|B_{\perp }\right|$ distribution will not deviate significantly from the designed configuration. This is verified by the calculated $\left|B_{\perp }\right|$ from simulated coil currents, as shown by the contour plots in Figure 9c–e: even at $t=16 \mathrm{m}\mathrm{s}$ when the total eddy current is~170 kA (see Figure 4d), the calculated field-null configuration still perfectly matches the designed one in Figure 6a, and the area of $\left|B_{\perp }\right|<$ 1 Gauss is even larger than that at $t=0 \mathrm{m}\mathrm{s}$. This demonstrates the high reliability of the field-null optimization method in Section 3.1, as well as the feasibility of the calculated coil current waveforms as control targets in real experiments. In addition, if the field-null configuration could be realized, it may potentially help mitigate the vertically unstable problem in early-stage plasmas due to the midplane-aligned breakdown region, which has a very low value of vertical elongation.

5. Discussion and Outlook

The simulation results shown in Figure 9 clearly demonstrate that stable and high-quality field-null configuration is not only designed but also has the potential to be realized in real experiments. This provides an exciting possibility of using this method on HL-3 to achieve more robust plasma breakdowns despite varying conditions, like pre-gas fueling, wall conditions and loop voltages. Yet, it is also clear that more verifications are required before actually using it in experiments with enough confidence. Firstly, the toroidally axisymmetric 2-D model of the vacuum vessel applied here could differ from the 3-D reality due to the existence of windows, pipes and first wall modules, etc. Thus, more careful calibration of the model parameters needs to be performed using supplementary diagnostic data, like magnetic probes and flux loops, etc. The construction of 3-D models, like using the finite-element method (FEM) [4], will be of great help in raising the accuracy in calculating passive eddy currents but is more complex than the presented model and is left for future work. On the other hand, the models of power supplies can also be improved to include the characteristics of thyristor rectifiers as they are on HL-3 [29], which could help in resolving the existing oscillation problems when performing FB control of CS and PF coils in experiments. Finally, targeted experimental tests are definitely needed to verify the optimized field-null configurations, as well as the applicability of designed coil current waveforms, and to further focus on varying operational scenarios and different boundary conditions to ensure that the configurations are robust and reliable under a wide range of conditions, if sufficient shots are available.

Since the $B_{\perp }$ configuration during the early stage not only affects the breakdown process but also impacts the equilibrium and position control during the subsequent plasma ramp-up stage, it is necessary to consider the optimization problem of $B_{\perp }$ in a more integrated way, as has been demonstrated in existing references. This could be performed in the proposed method by simply replacing the cost function $f\left(\mathit{I}\right)$ in Equation (8) with other forms. For example, the following $f\left(\mathit{I}\right)$ can be chosen to obtain a vertical field $B_Z$ closed to a wanted value $B_{Z,0}$ while retaining the field-null characteristic by minimizing $B_R$:

$$\begin{gathered} f\left(\mathit{I}\right)={\sum }_j^{j\in \mathbf{A}}\left|B_{R,j}^2\left(\mathit{I}\right)+{\left(B_{Z,j}\left(\mathit{I}\right)-B_{Z,0}\right)}^2\right| \\ ={\sum }_j^{j\in \mathbf{A}}\left\{{\left[{\sum }_{i=1}^{105}{G}_{R,ij}{I}_i\right]}^2+{\left[{\sum }_{i=1}^{105}{G}_{Z,ij}{I}_i-B_{Z,0}\right]}^2\right\} \end{gathered}$$

(17)

The optimization problem of $f\left(\mathit{I}\right)$ in Equation (17) is still quadratic and can be solved by transferring it into a linear equation set, similar to the process in Section 3.1. Other constraints, such as the magnitude of coil currents, could also be added in this way or by alternatively introducing boundaries to define the proper parametric domains. This indicates that the approach proposed in this paper is highly expandable and can be improved for various purposes in the future.

Real-time algorithms could also be helpful if more signals other than coil currents are involved, such as the data from magnetic probes and flux loops. By developing more complicated real-time algorithms, the control of plasma position or even shape could become possible in the early stages of plasma formation. Unfortunately, the circumstance here is complicated by the quick evolution of early-stage plasmas and large eddy currents in the vacuum vessel, and this part of work is left for the future.

Last but not least, we have noticed that artificial intelligence (AI) techniques based on various types of neural networks have been gradually applied to solving more and more real-time tokamak plasma control problems, including the prediction and mitigation of disruptions in advance [30,31], the accurate and flexible control of plasma boundary shapes [32] and sustaining high-performance plasma operations while getting rid of dangerous Magneto-hydrodynamic (MHD) instabilities [33], etc. Since the optimization problem of the whole tokamak plasma startup phase is generally integrated with some complex dynamics that are only exclusively known but hard to describe with high accuracy, it is actually similar to other control issues, like disruption prediction or control over MHD instabilities. Thus, it could be worth investigating the possibility of applying these powerful AI techniques to solve the tokamak startup design problem. One candidate here is the one utilizing the concept of “fuzzy” inference with the technique of fuzzy divergence, which exploits the ability of fuzzy logic to deal with uncertainty and non-linearities in complex systems. This technique has been successfully applied to both tokamak [34] and non-tokamak [35] problems, which shows fairly good transversality. Applying these kind of techniques to the problem in this paper could be significantly beneficial and is of great interest in future works.

6. Conclusions

In summary, this paper presents the work on building a complete set of simulation tools for optimizing the field-null configuration to achieve robust plasma breakdown on HL-3. First, a forward simulation model of the tokamak circuit is developed based upon the Simulink™ framework, including the coils, their power supplies and the passive conductive structures. The model has been calibrated, and the calculated coil currents can perfectly match the experimental values. With the developed model, an iterative approach is developed for optimizing the field-null configuration, which naturally includes the eddy currents in passive structures in a self-consistent way and can provide the CS and PF current waveforms that generate the optimum field-null configuration inside a given region. Finally, a PID feedback controller was designed and integrated with the model to form a closed-loop system, and simulation FB control experiments have demonstrated the feasibility of the designed field-null configuration and the corresponding coil current waveforms. These results build a solid base for applying the optimum field-null recipe to HL-3, which could help improve the quality and robustness of plasma breakdown on HL-3. This can help increase the machine productivity and could also assist in robustly achieving higher plasma currents by lowering the loop voltage at breakdown and saving flux consumption. The proposed method could also provide an alternative and more flexible way for designing the field-null configurations for future tokamaks like ITER.