{'fraction_non_alphanumeric': 0.024251189503028778, 'fraction_numerical': 0.009893311534511309, 'mean_word_length': 4.9222939424031775, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Development and Validation of a Tokamak Skin Effect Transformer model

A control oriented, lumped parameter model for the tokamak transformer including the slow flux penetration in the plasma (skin effect transformer model) is presented. The model does not require detailed or explicit information about plasma profiles or geometry. Instead, this information is lumped in system variables, parameters and inputs. The model has an exact mathematical structure built from energy and flux conservation theorems, predicting the evolution and non linear interaction of the plasma current and internal inductance as functions of the primary coil currents, plasma resistance, non-inductive current drive and the loop voltage at a specific location inside the plasma (equilibrium loop voltage). Loop voltage profile in the plasma is substituted by a three-point discretization, and ordinary differential equations are used to predict the equilibrium loop voltage as function of the boundary and resistive loop voltages. This provides a model for equilibrium loop voltage evolution, which is reminiscent of the skin effect. The order and parameters of this differential equation are determined empirically using system identification techniques. Fast plasma current modulation experiments with Random Binary Signals (RBS) have been conducted in the TCV tokamak to generate the required data for the analysis. Plasma current was modulated in Ohmic conditions between 200kA and 300kA with 30ms rise time, several times faster than its time constant L/R\approx200ms. The model explains the most salient features of the plasma current transients without requiring detailed or explicit information about resistivity profiles. This proves that lumped parameter modeling approach can be used to predict the time evolution of bulk plasma properties such as plasma inductance or current with reasonable accuracy; at least in Ohmic conditions without external heating and current drive sources.

treated using an additional non inductive current source [1], [2]. The coupling between transformer primary and plasma is accounted by a mutual inductance M. This standard description, despite its simplicity, has proved to be very useful in analysis and prediction of plasma discharge evolution. It has been used for years to infer non inductive current drive fractions and efficiencies in tokamaks [3], [4] and also to perform initial estimations of central solenoid magnetic flux in tokamak reactor studies [5] . However, it does not take into account the variations of plasma inductance that occur during plasma current transients due to the slow flux penetration inside the plasma, or skin effect. This is an important issue for tokamaks, particularly during the initiation and termination phases of the discharge, due to the large machine size and plasma temperature. We have developed an improved transformer model that includes a lumped parameter formulation for the skin effect, and whose purpose is to predict with reasonable accuracy the inductance and current changes as function of the external PF currents, plasma resistance and current drive. We call this model a skin effect transformer model to differentiate it from the standard transformer model where plasma inductance is a fixed parameter.
Distributed parameter simulations are the preferred option to simulate current profile evolution and the associated internal inductance changes that occur during plasma current transients [6] - [13]. These distributed parameter simulations predict strong correlations between internal inductance and plasma current transients, in agreement with the experiments. A correlation model is some kind of lumped parameter formulation, in the sense that does not contain explicit information about profiles. Since plasma current can be predicted with a lumped parameter model, the correlation observed between internal inductance and plasma current transients is evidence suggesting that a lumped parameter formulation could predict internal inductance behaviour. Previous lumped parameter modelling of the JET tokamak has proved the feasibility of this approach [14] .
To be able to describe the effect that plasma current profile evolution has on the internal inductance using the lumped parameter approach, some profile and geometric information must necessarily be embedded (lumped) in the model inputs, parameters and state variables. But the aim the lumped parameter modelling philosophy is not to describe plasma profiles or geometry evolution in itself. For this, distributed parameter models and simulations are required. There are plenty of such distributed parameter models available [15] - [21] , some of them are control oriented [22] , [23] and some of them are even available in real time [24] . From this point of view, lumped parameter models are no substitute for distributed parameter models. However, when it comes to control systems design, a lumped parameter formulation is preferred. Lumped parameter models have been the first modelling choice for tokamak magnetic control for years [25]- [30] , since control system's design methodology is well developed for lumped parameter systems [31] - [34]. The main motivation for the skin effect transformer model is to be used in the design of tokamak plasma current, boundary flux and inductance control systems. An ongoing project at TCV aims to test this type of controls. However Tokamak magnetic control is not the subject this paper. The reader is referred to [35] , [36] and references therein for details of how plasma current and general magnetic control systems are designed. The motivation, application and demonstration of inductance control are a reported elsewhere in the literature [37]- [42] and not discussed here either.
The paper is organised as follows: Section II starts the lumped parameter modelling by identifying some energy storing and dissipation elements of the tokamak circuit, along with some useful relationships among them. A set of exact coupled differential equations for internal inductance and plasma current are derived from energy conservation and power balance theorems. An additional boundary flux conservation theorem provides the link between external coil currents, plasma current and boundary flux. The result is a novel description of the tokamak transformer using only lumped parameters, predicting the evolution and non linear interaction of the plasma current and internal inductance as functions of the primary coil currents, plasma resistance, noninductive current drive and an additional equilibrium state derived from the flux linked by the plasma current distribution. This lumped parameter formulation does not require information about the plasma profiles spatial distribution. Instead, all the information about the profiles is embedded in the state variables and inputs.
Section III introduces a simplified model relating the equilibrium state to the boundary and resistive voltages. This is required to close the system of equations obtained in section II. The equilibrium state model is a lumped parameter version of the skin effect obtained by replacing the loop voltage profile eigenfunction by a three-dimensional eigenvector.
Section IV writes the equilibrium state model obtained in section III as a generic family of state space models, which is the preferred structure for the analysis of sections VII and VIII .
Section V explains the set-up and details of dedicated plasma current modulation experiments performed at the TCV tokamak. To expose the model non-linearities, the experiments were conducted using an extended range of operation both in amplitude and frequency. This differs from the conventional approach in system identification, where small signal models are sought. Some details of the TCV current control system are explained since they are central to the discussion, but this control system is part of earlier developments at TCV, and is not a contribution of this work.
Section VI explains how the data required for the analysis is generated with an equilibrium code, and some basic spectral analysis of the resulting inputs and outputs to give an idea of the experiment's modulation bandwidth.
Section VII performs a linear system identification analysis to find empirically the equilibrium state model order and parameters from the experimental data.
Section VIII performs a non linear system identification analysis using the model parameters found in section VII as an initial guess. A final cross-correlation study between the model residuals and its inputs is used to illustrate the accuracy of the results.
Finally, section IX presents a reduced version of the model written in control affine form, so it can readily be used for control systems design.

II. Derivation of exact tokamak transformer lumped parameter equations
We are after a lumped parameter formulation for the tokamak transformer not requiring explicit information about plasma profiles spatial distribution, but still being able to describe processes related to plasma profiles spatial evolution, such as the skin effect. Necessarily, some kind of profile information must be lumped in some states and inputs. This section identifies state variables and inputs that lump the profile information along with basic relationship among them. These relationships are derived from conservation laws, so the resulting model is exact.
In the following, both a cylindrical coordinate system ( ) , , r z φ and a flux coordinate system ρ,θ,φ ( ) are used, and the plasma is assumed to be axisymmetric about the z-axis. Only the time evolving components ( ) Z r B B , of the poloidal magnetic field and the toroidal components of the electric field E are considered in the analysis.
The region of integration will be defined as the region where there is plasma. This will correspond to a plasma volume G , or a plasma cross section Ω, delimited by the plasma boundary Γ . The plasma current is defined as where dS drdz = and j is the toroidal current density.
The poloidal flux function ψ r, z ( ) is the flux through an arbitrary circle of radius r centred at the torus symmetry axis at a position z. All the points with equal flux define a flux surface.
The flux surface surrounding the plasma region is the boundary flux surface ψ B .
The Green´s function for a unit toroidal current ring at r j = r j , z j The mutual inductances M j between the coil system j and the plasma boundary can be written in terms of Green's functions as [44] , [45] : Similarly, the external inductance between the plasma current and the plasma boundary is where dl θ ( ) denotes the differential arclength along the coordinate θ and Practical numerical details on how to calculate Green's functions and avoid the singularity at r = r 1 can be found in [44] , [45] Using the external and mutual inductances, the flux at the plasma boundary B ψ can be written as the sum of the contributions from the PF coil systems and the internal plasma current distribution, according to the boundary flux conservation equation: The plasma internal inductance i L is defined from the magnetic energy content inside the plasma volume G ( ) where 0 µ is the vacuum magnetic permeability and the differential volume element is The quantity W contains magnetic energy created by the plasma current as well as by external conductors.
The plasma inductance p L is the sum of the internal and external inductance components.
An alternative expression for the magnetic energy (7) as a function of poloidal fluxes at certain locations and the current density is [14] ( ) and equivalently ( ) where C ψ is defined as The flux C ψ is termed equilibrium flux, as it is the sole variable in our treatment that depends explicitly on the details of the plasma equilibrium. Although C ψ is a distribution average, we can formally identify a flux surface inside the plasma on which the flux takes the value C ψ .
This flux surface, which may of course move radially in time, will be termed the equilibrium flux surface.
The Poynting energy theorem in toroidal geometry is written as or equivalently as where the loop voltage at the boundary flux surface is (15) and the resistive drop R V is written in terms of the total plasma current, resistance and an equivalent non-inductive currentÎ as ( ) The resistive flux is defined as Again, we can formally identify a flux surface inside the plasma on which the flux takes the value ψ R . This flux surface, which may of course move radially in time, will be termed the resistive flux surface.
The ohmic power dissipation from Poynting analysis is given as a volume integral of the current density distribution j , effective plasma resistivity η , and non-inductive current density ĵ in the toroidal direction.
The ohmic power dissipation in the lumped parameter picture is R I −Î ( ) I = RI 2 − RIÎ (19) both descriptions are desired to be the same, so power balance is fulfilled: so the natural choices are leading to the following expressions for plasma resistance R and current drive Î : These are the only definitions that make compatible the Poynting power balance analysis (13) with (14) Linear combinations of (14) and of the time derivative of (11) lead to the following set of coupled equations ( ) The equilibrium voltage at the equilibrium flux surface is derived from Lenz´s law Equations (25) and (26) determine the evolution of the plasma current and internal inductance as functions of boundary size and shape, equilibrium, and resistive voltages.
The final steady-state solution for plasma current and inductance corresponds to a constant loop voltage profile across the plasma: The flux diffusion dynamics is embedded in this model. For instance, a drop in plasma resistance at constant boundary voltage will increase the plasma current, which in turn will increase the resistive drop (16). Similarly, an increase in boundary loop voltage will drive more plasma current and the resistive drop will also increase. Using the flux balance equation (6) and Ohm's law (16), the dynamic response (25), (26) can be written explicitly as a function of PF coil currents, plasma resistance and non inductive current drive: Combination of(6), (15), (29), (30) leads to an expression for the boundary voltage as direct function of the states, inputs and model parameters: The differential equations (29),(30) (or alternatively (25) , (26) , (31) ) predict the exact evolution of plasma current and inductance without requiring explicit information about plasma profiles spatial distribution. Instead, all the information about the profiles is lumped in the state variables, parameters and inputs. Note that (29), (30) is not yet a closed set of equations. Equations (29),(30) depend on the equilibrium state (27), whose dynamics will determined empirically in later sections using system identification techniques.
The plasma is kept within the vacuum vessel boundaries by means of additional plasma shape and position control systems acting on the PF coils. When shape and position are kept approximately constant, the terms on the right hand side containing time derivatives of mutual and external inductances can be neglected, and the set of equations simplifies to ( ) The physical meaning of this equation is that the change of flux produced by the external coils generates a loop voltage that compensates the resistive drop and builds up the plasma current. According to (34) , the time constant for inductive current build up is with R and L p given by (23) and (9).

III Lumped parameter model approximation for the skin effect
Equations (29),(30) are exact. They predict the evolution and non-linear interaction of the plasma current and internal inductance as functions of the primary coil currents, plasma resistance, non-inductive current drive and equilibrium loop voltage. To obtain a closed set of equations, we still have to obtain a model for the equilibrium voltage (27) state dynamics as function of the remaining system states and inputs of the model.
We want to find a lumped parameter version of the flux diffusion equation as function of the available states and inputs alone. The required output of the model is the equilibrium state. To build the lumped parameter equivalent to flux diffusion equation we have to determine what system states and inputs are relevant. To do this, we will replace the loop voltage profile eigenfunction by a three-dimensional eigenvector.
At any plasma location, the current density and vector potential are related by Poisson´s equation Writing the vector potential as a function of poloidal flux, Poisson´s equation can be written for the toroidal component of the current density as where * ∆ is the elliptic operator (2) .
Using Ohm's law, the voltage at any plasma location, and in particular C V , is a solution of Taking the time derivative of this expression leads to ( ) This distributed parameter equation is exact, but since the flux, resistivity and spatial distribution of non inductive current density are assumed to be unknown, it can not be directly incorporated in the model.
The simplest form of (39) corresponds to fixed plasma geometry with time invariant resistivity and non inductive current profiles; Thus, under these conditions, the time derivative of the loop voltage is proportional to the second spatial derivative of the loop voltage profile.
are the effective radii of the flux surfaces , , ψ ψ ψ , and the loop voltage profile varies smoothly across the plasma, we can approximate the second order spatial derivative by finite differences at neighbouring points around C ρ Equation (41) has already the shape of a stable first order differential equation for V C in which the resistive drop and the boundary voltage add up as a compound input to the system. It also has the required steady state solution (28). The metric (42) is unknown, and is also a function of time, as the equilibrium evolves. However, despite the numerous unknowns, (41) provides a useful insight into what states and inputs could be used in a linear state space model. This equation suggests in particular a state space model relating the equilibrium voltage C V to a linear combination of , We can also infer that, since the metric is expected to evolve in time, a first order structure with fixed coefficients is unlikely to be an accurate description, but we will not exclude this possibility from the start.

IV State space model for equilibrium voltage
For system identification purposes, we find convenient to write the dynamics of the equilibrium voltage (41) as a generic family of state space models of order n written in innovations form: Where x is a state space vector, is the input vector.The residuals r will represent all features present in the output that are not explained by the model. A steady state Kalman [46] matrix K is incorporated to reconstruct frequency response diagrams, including the noise power spectrum, perform residuals selfcorrelation (whiteness tests) and input-residuals correlations (independence tests). This will discussed in detail in later sections. The model order n and the matrices A, B, K can then be determined from experimental data using system identification techniques. For these techniques to work we need to dynamically stimulate the plasma to obtain loop voltage data with broad spectral content. To generate the required data for the analysis closed loop experiments using random binary signals (RBS) imposed as plasma current references have been conducted in the TCV tokamak. The preparation and details of these experiments are discussed in the next section

V System identification experiments
The 'Tokamak à Configuration Variable' TCV is a fusion plasma research facility particularly suited for the study of advanced plasma control systems [47]. It has a major radius of 0.9m, a typical plasma volume of ~1.5 m 3 , and a great flexibility of plasma shape configurations. TCV has an air core transformer with an available flux swing of 3. On TCV, the plasma current is made to follow a prescribed waveform by acting on the primary transformer coil voltage. Plasma current feedback is necessary to compensate for the changes in plasma resistance and non inductive current. For a detailed description of TCV hybrid control system the reader is referred to [48] Ideally, plasma current feedback control ought to be disabled to test the open loop response of the system. This is however not practical, since without current feedback the plasma current undergoes uncontrolled excursions, leading to plasma disruptions in many occasions. The situation is particularly acute in the extreme modulation conditions of the experiment discussed later in this section. We are then forced to maintain the plasma current feedback and adopt a closed loop system identification approach. Feedback systems introduce correlations between the perturbations in the feedback loop and the inputs. Unmodeled perturbations are treated as noise. When inputs and noise are comparable (as is the case when small signal models are being sought) closed loop system identification is better performed by adding a known perturbation inside the feedback loop [ [48] ], [49]. However, the high gain in the TCV current feedback loop implies that any perturbations introduced inside the loop (e.g. adding a feed-forward voltage at the ohmic coil, as in [49]) are quickly attenuated, so the perturbation technique was not practically compatible with the need to limit the current transients. Closed loop system identification, however, can be successfully performed when the system is stimulated with large signals compared with the noise input, and/or a good model for the noise input exists, so noise becomes a modeled perturbation. We are in fact interested in exploring the system's non-linearity in an extended range of operation both in amplitude and frequency, and the perturbation is largely coming from plasma resistance variations, for which we have a very accurate model at our disposal. This will allow us to use a direct closed loop system identification approach to extract the relevant open loop dynamic information, in both the frequency and time domains. To excite high order dynamics of the plasma current / inductance over a broad spectrum, RBS waveforms are imposed as references for the plasma current during the flat top. Plasma current was modulated between 200kA and 300kA with 30ms rise time, several times faster than its time constant L p /R≈200ms. A feedback control system generates primary transformer voltage control actions to drive the desired plasma current. The plasma current feedback system has a settling time of less that 8ms, negligible compared with the time constant (35) τ p ≈ 200ms . The resulting voltage signals for the primary transformer coil and the corresponding responses in primary coil current and plasma current are shown in figures 1 and 2. The time derivative of the external flux created by the transformer coil and the additional PF coil currents create a plasma voltage that for constant plasma and mutual inductances can be approximated by (34). The primary transformer current excitation and the resulting plasma current are shown in Fig. 1 for two system identification experiments. A positive current in the toroidal direction flows in the direction of increasing toroidal angle. The sign criteria are given by Lenz's and Ohm's law. A positive voltage applied to the transformer primary generates a current that increases in time during the discharge. The corresponding boundary flux that increases in time will generate negative boundary voltage and negative plasma current. The time prior to the plasma current build up (shown as negative) is used to charge the primary transformer with negative current, to increase the available flux swing. Then, the transformer current is brought up to positive values to drive negative plasma current. The time origin 0 0 t = is chosen as the time when the poloidal magnetic field null necessary for breakdown is created inside the vacuum vessel. Figure 2 shows the same quantities during the RBS time window, with the offsets and linear trends removed.

VI Preliminary data generation and analysis
The required inputs and outputs to the system cannot be measured directly. The data used for the validation is obtained from the equilibrium reconstruction code LIUQE [50] . This code determines magnetic flux surface distribution from a set of external magnetic measurements [51] . Among other plasma variables, LIUQE delivers the plasma current, boundary flux and normalised internal inductance that are required for the analysis. Plasma current and internal inductance can also be obtained from magnetic probes surrounding the plasma [52] but we prefer to use the LIUQE's plasma current to obtain a consistent data set. Similar results can also be obtained using transport simulations [53] .The plasma internal inductance is defined relative to the machine major radius ( 0 0.88 r m = on TCV) and the normalised internal inductance, using The resistive voltage is derived from the power balance (14) as The equilibrium flux is derived from (11) as Boundary and equilibrium voltages are obtained from (15) and (27) respectively.
External magnetic measurements from flux loops and probes around the vacuum vessel are acquired at a 2 kHz sampling rate. The reconstructed equilibrium data for the current analysis has been generated with a 1.5 ms spacing.  figure 5 . The spectral power density of the signals halves above 30 Hz respect to its maximum. Above 70 Hz there is no significant power in inputs or outputs, so it will be very difficult to reconstruct the spectrum (bode diagrams) above this frequency. This will be shown in the next section.  Power spectra of input u=0.5(V B +V R ) and output y=V C signals for two system identification experiments. Both system identification experiments are almost identical in terms of power spectral content. Most of the power is concentrated below 70 Hz.

VII Linear system identification
Next, we conduct a linear system identification analysis from actual experimental data, to determine the best order of the model. To do this, the experimental data from a RBS stimulated discharge (39950) is used to fit the model parameters using the prediction error method. The parameters found are then used to predict the behavior of another discharge (39952) with an infinite prediction horizon (simulation). Finally, the equilibrium voltage simulation is compared with the actual experiment, and a loss function R N is derived from the residual ε .
where N is the number of experimental y x and simulated y s data points in the prediction set. The sample variance of the data set is defined as and the overall fit parameter is defined as The Akaike Information criteria (AIC) and forward prediction error (FPE) are used to measure the quality of the prediction and best system order [54]. AIC and FPE not only reward the goodness of fit, but also include a penalty that is an increasing function of the number of estimated parameters. According to Akaike's theory, the most accurate model has the smallest AIC and FPE. The AIC and FPE are defined by the following equations: where d is the number of estimated parameters, and C is a constant that that be ignored in model comparisons [54]. Figure 7 shows the simulation results for models up to seventh order. The fit parameter for the first order model (red trace) is 1.5%. Clearly, a first order system (red trace) is not sufficient to describe the dynamics. This supports our initial insight, stemming from the observation that the metric of (41) evolves in time along with the loop voltage profile. Generally speaking, the higher the order the lower the loss function will be, at the expense of introducing a larger number of parameters in the model. Despite the good visual agreement of most models of fig. 6, the fit parameter (51) for the second order system simulation is just 25%. This is explained by the small scale noise content features in the identification data set, which are difficult to appreciate in the figure. Next we evaluate the impulse response of the system using a non-parametric model [54] derived from the experimental data using correlation analysis. Because nonparametric and parametric models are derived using different algorithms, agreement between these models increases confidence in the parametric model results. The corresponding transient plots provide also valuable insight into the characteristics of model dynamics, such as natural resonances or settling time. We consider the simplest, discrete-time model to estimate the dependence of current and past output values subject to noise in the system. We estimate a high order (n=1 , m=70 ) autoregressive exogenous (ARX) model [54] , [55], [56] with the structure:

y t a y t T a y t n T b y t b y t T b y t m T e t
A small negative lag is introduced in the correlation to investigate any residual feedback effects in the input-output data. Then we use this high order model to simulate the system's impulse response. Finally, the parametric models up to seventh order are compared with the correlation model response. Figures 8 and 9 show the results. The first order model is the worst performer. The second order approximation seems closer to the under-damped behavior and natural frequency (about 30Hz) shown by the correlation model.
Good models have the self correlation of the residuals (whiteness test) and the cross correlation of the residuals with the input (independence tests) below certain confidence intervals that are calculated from the estimated uncertainty in the model parameters, assuming a Gaussian distribution for the estimates. The results of this analysis are shown in Fig. 10 for models up to order 3. The correlation windows are extended over one characteristic response time of the system (35), the order of 200 ms. The first order system is the worst performer in this set. Second, third and higher orders (not shown) pass the whiteness test, but only the third order passes the correlation test, followed marginally by the second order model. Fig. 10. Whiteness and independence tests for models up to third order. Self correlation of residuals (top) and cross correlation between input and residuals (bottom) are shown for up to 75 lags (112.5 ms). The discontinuous lines bound the range that has a 99% probability of being statistically insignificant. Figure 11 shows a comparison between the model and experimental data frequency responses. The experimental frequency response is not well determined above 50Hz, particularly in phase, due to the low signal to noise ratio at these frequencies. This can be appreciated in figure 5, where it can be observed than the spectral power density of the input signals halves above 30 Hz respect to its maximum. The first order system response is completely inadequate, whereas the second order response seems again to be closer to experimental results, exhibiting a resonance frequency around 30 Hz. The second order system with two poles and one zero is consistent with the experimental data, despite the noise present. The 45°/decade rise of the zero is arrested by the overlapping of the two poles, each contributing with a 45°/decade drop in the limited range of frequencies where all are active contributors to the phase. According to all the tests, a second order linear model seems to be the best approximation for equilibrium loop voltage dynamics. The resulting second order state space model can be described by the following transfer function: The parameters found for the second order system using the previous analysis are summarised in table 1 ω δ T k 174 0.3 0.012 0.9 Table 1 Linear system identification of a linear second order model using the equilibrium voltage V C in the objective cost function.
We generally trust the model order found by the system identification, and believe the frequency ω at which the equilibrium loop voltage will resonate when driven by an external oscillation (see table 1) is probably in the right magnitude order (tenths of Hz). The steady state gain found (k=0.9) is also quite close to what we expect in steady state (k=1). But we don't trust the damping factor found, for the following physical reasons. The damping factor is 0.3 δ = , allowing resonant natural oscillations on the equilibrium loop voltage even when the system is not driven by an external source. This resonant behaviour requires at least two forms of energy accumulation in the system that can mutually be exchanged. In an electric circuit these are the energy stored in the electric fields (e.g capacitors) , and the energy stored in magnetic field (e.g. inductors). In our problem the energy stored in the poloidal magnetic field is about 20 KJ, and the energy stored in the toroidal electric field (<1V/m) is negligible compared with it (<10 -11 J) . The time variations of toroidal electric field are associated to the displacement current. At frequencies below the GHz range these can be neglected in the tokamak analysis. Therefore, if the toroidal electric field was responsible for the effect, resonant behaviour should be observed at frequencies above the GHz range, not in the sub-kHz range as we observe. Another possibility is the coupling with the radial electric fields through the return current [57]. In this case the radial electric fields are the order of kV, but again the energy contribution is again too small (<10 -8 J) compared with the magnetic energy. Yet one more form of energy accumulation is kinetic in nature. During the experiments, plasma radial excursions take place as a response to the internal current redistribution. In this experiments this can take place at maximum speeds of 0.5 m/s. This gives a plasma kinetic energy term [28] which is again negligible (<10 -7 J) compared with the energy stored in the magnetic field, so it could not be responsible for the resonance found. None of the above forms of secondary energy accumulation can explain the resonant behaviour in the sub-kHz range. We should have obtained instead a damping factor larger than one, corresponding with a damped system δ ≥ 1 with two real poles (two time scale dynamics), which does not allow oscillations unless the system driven by an external source.
We believe the observed resonant behaviour may be due to artefacts introduced by the equilibrium reconstruction, or perhaps due to the fact that we have used a single compound input for the analysis. There are also inputs to the system that have been assumed to be small and have not been modelled in the equilibrium voltage dynamics either, like bootstrap current. There are also conceivable uncertainties in the equilibrium reconstruction that are not considered, besides of course, all features that can not be explained due to the limitations of the lumped parameter approximation. All these factors are considered to be small contributors, and are included in the model as small perturbations or noise, which are accounted by means of the Kalman noise model structure of (43). The noise spectrum of these is presented in figure 12 exhibiting a maximum that suspiciously coincides with the found system resonance (27Hz). Compared with its maximum values (see figure 5), the input and output signals at the resonance have about four times less power, translating in about 16 times less signal amplitude. So another possibility is that the relative contribution from the true system dynamics and noise are difficult to resolve by the linear system identification, due to the small signal to noise ratios at this frequency range. To obtain more accurate parameters we will have to resort no non-linear system identification techniques.

VIII Non linear system identification
In this section we are going to exploit the full non-linear model to identify the parameters of the equilibrium voltage model more accurately. We will build a cost function based on the triplet , , i C L I V and exploit the constraints imposed by the model's non-linear dependencies. We start by using separate inputs , B R V V and restrict the state space model to a second order structure without steady state Kalman gain. To be able to determine the most appropriate set of initial conditions, the second order approximation for equilibrium loop voltage dynamics that was found in the previous section is written in canonical observable form as ( ) ( ) and ω ,δ are the natural frequency in radians per second δ damping factor. The variables have been scaled to the natural frequency so that the model parameters are all roughly of the same magnitude. Written in this way, the corresponding transfer functions are The difference with the model found in the previous section is that now we will allow separate inputs for boundary and resistive voltages. Choosing a relative weight 1 k ≤ for B V , and a relative weight ( ) 1 k − for the R V input, the steady state solution of (58) is automatically guaranteed to match the steady state solution (28).
The weight k regulates the relative contributions from boundary and resistive voltages to the equilibrium voltage. For instance, if k = 0.5 both inputs account for an equal 50% contribution. This weight should be approximately in inverse relationship with the relative proximity of the equilibrium and boundary flux surfaces, according to (41).
The natural frequency ω and damping factor δ are a generalization of a "two time scale" model, which allows for resonant behavior when the damping factor is less than one. This allows us to explore new physics such as the existence of a plasma capacitance that could give resonant behavior on the loop voltage at high frequencies. As it was discussed in the previous section, this type of resonant behavior should not be observed at the sub-kHz range, but the linear system identification of the previous section has given some importance to this possibility, so we should retain it. When the damping factor is greater than one, the system has two real poles, and a typical "two time scale dynamics" is recovered. Regardless of the damping factor value, the natural frequency should scale in inverse relationship with the plasma temperature, just like the skin time scales in direct relationship with temperature. The factors T B , T R , account for fast events in the high frequency end of the spectrum, so the equilibrium loop voltage knows about high frequency changes at the boundary or resistive voltages at faster rates than the resistive time scale. This accounts for processes like fast magnetic reconnection, turbulence, helicity transport [58] etc To find the model parameters, the linear diffusion model (58) is combined with the non linear equations (16), (25)and (26) to form a closed set of equations. The input R V to the linear diffusion model (58) has now become a non linear function of the plasma current, plasma resistance and current drive, according to (16). Equations (16), (25), (26) , (58) are then evolved in time starting from given initial conditions. The simulation results for both experiments (39950 and 39952) are finally compared with the observed data and a cost function is built as a weighted sum of the least squares prediction errors. A gradient descent algorithm is used to find the model parameters for which the cost function has a minimum.
Since the system has outputs measured in many different units at different magnitudes, volt, µH and MA, the weighting function is chosen as the inverse of the data sample variance (50) to regularize the units and avoid excessive weighting on some of the outputs. As with any nonlinear optimization algorithm, there is a chance that the model might find a local minimum that is not accurate. To help the non linear optimization process, an initial guess for the model parameters obtained from the linear system identification analysis presented earlier is used.
The results from this non-linear optimization are shown in figure 13 (red traces), and the corresponding parameters are summarised in table 2. Shots The parameters found using the non-linear optimization correspond again with an underdamped model. We would like to exclude this result because, as mentioned earlier, resonant behaviour cannot physically exist in the sub-kHz range. So we perform a final check restricting the damping factor to be 1 δ ≥ in the optimization process. Both shots are identified separately both to increase the agreement with experimental data and to check parameter value dispersion. The simulation results with the optimized parameters (blue traces) are shown in figure 13 The parameters found for the two shots along with their standard deviation are given in table 3. Shot (3) 1 (0.04) 0.027 (0.001) 0 (1e-5) 0.6 (0.01) 39952 440 (4) 1 (0.05) 0.028 (0.001) 0 (1e-5) 0.5 (0.01) Table 3 Non linear system identification using linear (V C ) and non linear (I, L i ) outputs in the objective cost function. The parameters are found for a damped second order model .
Looking at the results of figure 13, it seems the damped models of Table 3 (blue traces) are also compatible with the experimental results (black traces), with small differences with respect the under damped model of Table 2 (red traces). The overall fit parameters for plasma inductance (38%), current (70%) and equilibrium voltage (30%) are almost identical for both the damped and the under damped models. This evidences a flat cost function vs. the model parameters, with broad minima. We may inquire whether the low fit parameters obtained, particularly in the inductance and equilibrium voltage, are the result of small scale noise or represent a more fundamental failure to describe the dynamics. To this end, we will perform a correlation analysis between the model residuals and inputs. If the unexplained residuals correlate with the inputs, there are features in the outputs that can still be related to the inputs, and therefore these are unexplained features not reproduced by the model. If the correlation levels are below the confidence intervals, the discrepancy can be safely attributed to noise, and this would mean that we have reached the best possible model that can be extracted from the given data. The results can be found in figure 14. The cross correlation between the inputs and residuals for equilibrium voltage lay outside the 99% confidence intervals (dotted lines). There are features in the equilibrium voltage that correlate with inputs, and therefore are not explained by the equilibrium voltage model. The left-hand column shows the correlation functions between the plasma resistance input and the residuals from plasma inductance, current and equilibrium voltage. The right-hand column shows the correlation functions between boundary loop voltage and the residuals from plasma inductance, current and equilibrium voltage. The analysis is performed up to ±75 lags (±112ms). The discontinuous lines bound the range that has a 99% probability of being statistically insignificant.
We have no uncertainty information in the equilibrium reconstructed data, so we can not be sure about the accuracy of the second order model parameters found. Quite likely different equilibrium reconstruction methods could give different order and model parameters. Nevertheless, the lumped parameter model presented explains the most salient features of the plasma current transients, such as the inverse correlation between plasma current ramp rates and internal inductance changes, without explicit or detailed knowledge about current distribution, resistivity or loop voltage profiles. This proves that explicit or detailed knowledge of plasma profiles is not necessarily required to predict the time evolution of bulk plasma properties such as plasma inductance or current with reasonable accuracy; at least in Ohmic conditions without external heating and current drive sources. The reason for this is the smoothness of the loop voltage profile under these conditions. The skin effect transformer model provides a good balance between model accuracy and simplicity in this case.

IX Reduced model in control affine form
We consider the case of interest where plasma resistance and non inductive current are model parameters or perturbations, and the only actuators available are the currents in the OH coil and PF systems. We define the system input, output and state vectors as ( )   can be obtained as a combination of the system states as Written in the control affine form, the model is ready to be used for the design of control systems for plasma current, inductance or any other variable of interest that can be obtained as function of the state vector (61).

X Conclusions
A lumped parameter, state space model for the tokamak transformer including the slow flux penetration in the plasma (skin effect transformer model) has been developed and validated with TCV tokamak discharges.
The model does not require detailed or explicit information about plasma profiles or geometry. Instead, this information is lumped in system variables, parameters and inputs. Exact expressions for these and fundamental relationships among them have been derived from basic electromagnetic theory.
The model has an exact mathematical structure built from energy and flux conservation theorems, predicting the evolution and non linear interaction of the plasma current and internal inductance as functions of the primary coil currents, plasma resistance, non-inductive current drive and the loop voltage at a specific location inside the plasma termed the equilibrium loop voltage.
Loop voltage profile in the plasma is substituted by a three-point discretization, and ordinary differential equations are used to predict the equilibrium loop voltage as function of the boundary and resistive loop voltages. This provides a simplified model for equilibrium loop voltage evolution, which is reminiscent of the skin effect. A systematic procedure has been applied to find the order this differential equation and its parameters in terms of objective cost functions, quantitative information criteria, frequency response, and residuals whiteness/correlation tests.
Fast plasma current modulation experiments with Random Binary Signals (RBS) have been conducted in the TCV tokamak to generate the required data for the analysis. Plasma current was modulated in Ohmic conditions between 200kA and 300kA with 30ms rise time, several times faster than its time constant L p /R≈200ms. A second order differential equation has been found to be the best approximation for equilibrium loop voltage model. When the equilibrium voltage model is combined with the energy and flux preserving equations, the resulting non linear model is capable to describe the plasma current and internal inductance time evolution with 70% and 38% accuracy under these modulation conditions. The skin effect transformer model explains the most salient features of the plasma current transients, such as the inverse correlation between plasma current ramp rates and internal inductance changes. This proves that explicit or detailed knowledge of plasma profiles is not necessarily required to predict the time evolution of bulk plasma properties such as plasma inductance or current with reasonable accuracy; at least in Ohmic conditions without external heating and current drive sources. The reason for this is the smoothness of the loop voltage profile under these conditions. The skin effect transformer model provides a good balance between model accuracy and simplicity in this case.