Application of the Waveform Relaxation Technique to the CoSimulation of Power Converter Controller and Electrical Circuit Models
Abstract
In this paper we present the cosimulation of a PID class power converter controller and an electrical circuit by means of the waveform relaxation technique. The simulation of the controller model is characterized by a fixedtime stepping scheme reflecting its digital implementation, whereas a circuit simulation usually employs an adaptive time stepping scheme in order to account for a wide range of time constants within the circuit model. In order to maintain the characteristic of both models as well as to facilitate model replacement, we treat them separately by means of input/output relations and propose an application of a waveform relaxation algorithm. Furthermore, the maximum and minimum number of iterations of the proposed algorithm are mathematically analyzed. The concept of controller/circuit coupling is illustrated by an example of the cosimulation of a PI power converter controller and a model of the main dipole circuit of the Large Hadron Collider.
1 Introduction
Superconducting magnets are used in modern highenergy particle accelerators to steer the trajectory of beams of charged particles. The current in the magnets has to follow the increase of the energy of particles. A dedicated closed feedback loop is implemented for the magnets to generate the desired magnetic field. On top of that, there are quench protection systems responsible for discharging the stored magnetic energy in the circuit in case of emergency. A quench is a local transition from a superconducting to a normal conducting state. After a quench, the magnetic energy is dissipated in a small fraction of the coil and without any countermeasures, such a fault would lead to irreversible damage of the magnet and consequently the circuit. Computer simulations play an important role in understanding these complex phenomena [1]. Recently, a waveform relaxation algorithm has been implemented in the STEAM (Simulation of Transient Effects in Accelerator Magnets) project to perform field/circuit coupling for the magnetothermal analysis of superconducting magnets and circuits [2]. The goal of the present work on waveform relaxation for controllercircuit coupling is to include controller models in the STEAM framework. This will allow us to study the circuit behavior on a wider range of potential failure scenarios.
In this paper we focus on the nominal operation of a circuit, coupling models of a power converter controller and circuit. We consider a PID type current controller of a power converter and a linear electrical circuit. Typically, the respective models of both, controller and circuit, are developed separately with different software packages. There are several motivations for that. Firstly, both models are solved with different timestepping, i.e., the controller model is executed with fixed time step, whereas the circuit model employs an adaptive timestepping algorithm in order to account for the transient behavior. Secondly, the controller is usually designed on the basis of the firstorder equivalent electrical circuit models, similarly, a circuit driven by an ideal current source is satisfactory for a wide range of simulation scenarios. Nevertheless, a coupling between both models can provide more information on both controller and circuit behavior, especially for the relevant (nonlinear) failure scenarios, where ultimately partial differential equations are used to mathematically describe the quench.
In [3] authors employ waveform relaxation method to perform controller testing by means of simulations. This approach is an alternative to hardware in the loop simulations, however does not discuss the convergence properties. We propose the application of waveform relaxation [4, 5, 6] to the coupled controller/circuit problem and study its features. Instead of sequentially executing controller and circuit simulator according to the sampling frequency, the problem is translated into an iterative process on larger time windows. Such an approach promises benefits in terms communication costs and thus eventually computational time. Ultimately the approach will allow to seamlessly include the action of a controller in a STEAM fieldcircuit simulation; see above.
The remainder of the paper is organized as follows. Section II introduces the circuit and controller models. Section III discusses the proposed waveform relaxation algorithm and studies convergence properties. In Section IV we present a PI controller design procedure for the firstorder model of the Large Hadron Collider (LHC) main dipole circuit. The paper is concluded with an application of the waveform relaxation scheme to the simulation of electrical transients in the superconducting circuits powered by a voltage source controlled by a PI controller.
2 Mathematical Models
In this section we present a general setting for treating controllers of PID class and linear circuits composed of resistors, inductors, and capacitors as ordinary differential equations (ODEs) and differentialalgebraic equations (DAEs) on a time interval .
2.1 PID Controller Model
The governing equation of a continuous PID controller, with , , and denoting, respectively, proportional, integral, and differential gains, is given in continuous form as
(1) 
where denotes a possible delay due to physical limitations of the controller, and is the error given as a difference between the reference signal and its actual waveform . The temporal discretization is given by times with . The rectangle method for quadrature and the ImplicitEuler scheme for differentiation yield
(2)  
For the remainder of the paper let us consider the special case of a PI controller in order to analyze the convergence properties of the considered waveform relaxation method. Disregarding differential gains () and differentiating equation (1) yields the following (delay) differential equation
(3) 
which is trivial in the sense that the right hand side does not depend on the unknown. It can be rewritten as
(4) 
with the initial value , the unknown , and the excitation .
2.2 Electrical Circuit Model
Consider an electrical network composed of nodes with denoting the th nodal voltage with respect to ground, and branches with and being the current through the th branch and the voltage across the th branch, respectively. The circuit topology is described by the incidence matrix , with nonzero elements equal to either 1 or 1 depending on the assumed direction of branch with respect to node .
The incidence matrix allows writing the Kirchhoff Current Law as , where is a vector of branch currents. Furthermore, the vector of nodal voltages is related to the vector of branch voltages by .
In case of the LHC main dipole circuit during nominal operation, i.e., when the power converter is active, only linear capacitors, inductors and resistors are considered, and the incidence matrix A can be decomposed into block form
(5) 
where is the resistance incidence matrix, the capacitance incidence matrix, the inductance incidence matrix, and is the PID controller (voltage source) incidence matrix. Based on those matrices, a socalled conventional formulation of the Modified Nodal Analysis (MNA) [7] is expressed as
(6a)  
(6b)  
(6c) 
with capacitance, inductance, and conductance matrices , , and , respectively. The unknowns of the system are the vector of node potentials , and currents through inductors . In order to obtain a unique solution, an additional equation for the PI controller’s unknown current in terms of the applied voltage drop is added and a grounding node is selected.
In analogy to (4), the system (6) is reformulated into an abstract differentialalgebraic initialvalue problem as
(7) 
with initial value and unknowns . The external input function is given in the studied case by the output voltage of the controller which is obtained from ; see (4).
Let us discuss the simple case of an inductance connected in series to a resistance and a controller (see Figure 1); the incidence matrices read
with the elements and . The MNA yields four DAEs in the (redundant) variables , , and . It is mathematically equivalently expressed by the ODE
(8) 
Recovering the underlying ODE is straightforward for small systems and if there are no LI cutsets or CV loops present in the circuit [8].
3 Analysis of the Waveform Relaxation Scheme
In practice, the circuit model (7) is usually numerically solved with an adaptive time stepping algorithm (e.g. trapezoidal rule or BDF schemes, [9]) whereas the controller model employs fixed time steps corresponding to its sampling frequency. Coupling of both models via waveform relaxation allows to maintain this setup, by exchanging waveforms at communication time points . For that purpose, the overall simulation time is divided into time windows with . For each time window the two models are solved separately and their solutions (waveforms) are exchanged. The exchange of the waveforms can be organized such that (i) both models are provided with results from the previous iteration, (ii) one model is provided with a solution from the previous iteration and the other one with a solution calculated in the current iteration. The first method, called Jacobi scheme, is profitable for a parallel execution of models, whereas the second one, the GaussSeidel scheme, brings the benefit of faster convergence but it implies sequential execution of models [6]. The termination criterion for the waveforms’ exchange, i.e., for the convergence of the iterations in , is determined by applying an appropriate norm to measure the difference of two subsequent waveforms obtained with one or both of the models.
We employ the GaussSeidel scheme for the waveform relaxation method. At each time window, controller and circuit models are solved sequentially as depicted in Fig. 2. In order to incorporate the behavior of the delay a zero order hold (ZOH) system is present in the main control loop [10].
In the case of the considered coupled problem, with and denoting the time window and convergence iteration index, respectively, the waveform relaxation algorithm takes the following steps:

Set and extrapolate and for .

Solve (4) with initial value and input using fixed time steps to obtain for .

Solve (7) with initial value and input using an adaptive time stepping scheme to obtain for .

If not converged, then set and go to Step 1). Otherwise set and
If proceed Step 0) to start the next time window.
In the important case of ODEs, i.e., if there are no algebraic constraints
with , being piecewise continuous functions and being globally Lipschitz continuous with respect to for all , the GaussSeidel waveform relaxation algorithm converges uniformly on bounded time intervals (see [4, 6]). This implies that the controllercircuit waveform relaxation scheme converges provided the circuit is given as the ODE (8). This analysis can be extended to DAEs, e.g. [11, 12] and even to delay differential equations [10]. However, the classical reasoning disregards errors due to time stepping since this error can be made arbitrarily small by reducing the step size in the respective models. This is not possible in the case of an external controller, where the time step is not given by numerical concerns but dictated by the actual hardware. The following analysis discusses this timediscrete setting.
Proposition 1.
Proof.
Being and at , a fine enough discretization of (8) is assumed such that the time stepping error can be disregarded which yields
(10) 
with operator . For equation (3) an implicit Euler scheme is applied at time steps . This, combined with equation (10) leads to
with . It is sufficient to show that for , , as then , which implies the proposition. This follows via induction. ∎
An additional termination criterion can be introduced. The iteration stops when either (i) a difference between two subsequent waveforms is below a certain tolerance, or (ii) the convergence iteration index reached . Unless the current converged in less than iterations, we can observe that the iterationbased stopping criterion requires one convergence iteration less with respect to the currentbased convergence criterion. Additionally, if the time window size is equal to the controller sampling period, , according to Proposition 1, the algorithm reduces to a weak coupling scheme, i.e. both models are executed only once () and the convergence check is not carried out.
For a certain class of open loop systems and reference profiles, the number of convergence iterations in the steady state is less than as shown in Proposition 2.
Proposition 2.
For systems of type p, i.e. systems with pfold pole equal to zero in the continuous open loop system transfer function, with the reference signal being a polynomial function of time of (p1)th order, and for a given time window in the steady state, the proposed waveform relaxation algorithm converges at .
Proof.
The steady state error for systems of type p with reference signals that are polynomials of time of order p1, is zero
(11) 
where and are continuous transfer functions of explicit, input/output forms of (4) and (7), respectively. Since the steady state error is equal to zero, according to (2), the controller output does not change . The extrapolated controller output results in identical subsequent current profiles , which concludes the proof. ∎
Proposition 2 indicates the minimum number of convergence iterations, and together with Proposition 1 provides upper and lower bounds for the number of convergence iterations during a transient. We obtain .
The proposed waveform relaxation algorithm is now applied to a practical example of the cosimulation of LHC main dipole circuit and its powerconverter controller. The next section deals with a selection of a PI controller gains based on an equivalent firstorder approximation of the RB circuit.
4 PI Controller Design for a First Order Model of the LHC Main Dipole Circuit
The goal of the current controller is to follow a predefined current profile corresponding to the energy increase of the particles in the accelerator. In other words, the increase of the magnetic field in the main dipole chain is synchronized with the increase of the energy of particles provided by the accelerating cavities. To ensure continuity of the first order derivative, i.e., , parabolic joints are introduced between linear sections whenever the current derivative is changing its value.
For the controller design, we consider the LHC main dipole circuit, composed of 154 dipole magnets connected in series along with protection and powering devices. The main dipole circuit is represented as a first order system with equivalent inductance and equivalent resistance . The design will be carried out in frequency domain by means of the Laplace transform for the closed loop system shown in Fig. 3.
The transfer function of the firstorder circuit model is given as
(12) 
and the PI current controller transfer function reads
(13) 
The open loop transfer function is
(14) 
The characteristic polynomial of the closed loop system has the following form
(15) 
The selection of the controller influences the roots of the polynomial and in consequence the dynamic response of the system. For that purpose we compare (15) with the normalized second order polynomial equation
(16) 
with representing the relative damping, and where is the undamped natural frequency of the system. By comparing the coefficients of the respective polynomials (15) and (16) we derive the relations for the controller gains
(17a)  
(17b) 
For a practical application, the circuit should smoothly respond to a change in the reference current and oscillations are not desired. The typical value of relative damping is . The value of natural frequency is determined by the bandwidth of the reference current profile and is equal to leading to .
5 Numerical Examples
The performance of the proposed waveform relaxation scheme is verified by means of two simulation scenarios compared to a monolithic reference simulation. For circuit and controller simulation we employ ORCAD Cadence PSpice and Powersim PSIM, respectively. The monolithic simulations are carried out with Powersim PSIM. According to the procedure described in Section IV (17), the controller gains are equal to and . The controller sampling period is is chosen as controller step size and the maximum time step size for the adaptive timestepping algorithm of the circuit simulation is set to with absolute error tolerance equal to .
5.1 Step Response of a FirstOrder Model
In the first test we verify the algorithm operation and validate Proposition 1 and 2. We consider a step response of the firstorder model of the RB circuit cosimulated for s. The time interval is divided into windows of fixed length (). Communication between the solvers occurs at discrete time instants . In order to demonstrate Proposition 1 and 2 we consider only the currentbased termination criterion given as
Figure 4 shows the comparison between the reference step current profile and the current responses of the cosimulated and monolithic simulations. The test shows a good agreement between results obtained in both cases.
Controller output voltages for the first time window are given in Table 1, corresponding to iteration steps, as well as to the monolithic simulation. As one can notice, at each consecutive iteration, the controller output approaches the monolithic solution . In other words, each iteration extends the interval of time for which the controller output is equal to the monolithic reference solution, which is the typical behavior of the waveform relaxation scheme [5]. It is also worth noticing that the controller output for the last two iterations is identical, as the currentbased convergence criterion requires two subsequent iterations with the same voltage inputs.
The number of iterations for each time window are reported in Fig. 5. During the initial transient, the number of convergence iterations is equal to , whereupon it gradually decreases to two iterations, as the integral part of the PI controller calculates appropriate output for the studied circuit topology and excitation function. (Note that is the first iterate.) In this case, the open loop system is of type 1 as one can notice from (14). To conclude, the test underlines the correctness of the number of iterations predicted by Propositions 1 and 2.
[s]  0.04  0.08  0.12  0.16 

[V]  161.16  119.34  76.12  41.67 
[V]  161.16  185.48  209.80  234.11 
[V]  161.16  119.34  62.55  14.95 
[V]  161.16  119.34  76.12  44.45 
[V]  161.16  119.34  76.12  41.67 
[V]  161.16  119.34  76.12  41.67 
5.2 ParabolicLinear Response of a Realistic Model
After verifying the operation and convergence properties of the developed waveform relaxation coupling scheme, in the second test we perform a cosimulation of the considerably more complex LHC main dipole circuit composed of several thousand lumped components [13]. The underlying ODE of this realistic model cannot be easily extracted and therefore this setup may not be covered by classical arguments of waveform relaxation analysis.
The circuit is composed of 154 equivalent RLC models of a dipole magnet along with protection devices, a power converter and a filter. Parameters of the equivalent models were identified to match the frequency behavior of each of the magnets [14]. We consider nominal operation of the circuit, i.e., quench protection systems are deactivated. In this test, we verify the algorithm’s operation as a weak coupling scheme. The cosimulation is performed for s and covers the initial parabolic part of the current profile, concatenated with the linear increase. The time interval is divided into windows of fixed length (). The calculated current response accurately reproduces the reference current profile (see Fig. 6). The power converter controller output voltage is shown in Fig. 7.
6 Conclusion
In this paper, a waveform relaxation algorithm for the cosimulation of powerconverter controller and electricalcircuit models has been proposed. The algorithm divides the total time interval into windows, and for each window both models are solved separately with an appropriate timestepping algorithm. Properties of the waveformrelaxation algorithm’s convergence have been studied and stated in two propositions. The first proposition defines the maximum number of convergence iterations leading to an additional termination criterion and potential reduction of the number of iterations with respect to the currentbased stopping condition. The second proposition analyses the conditions for the minimum number of iterations. Boundedness of the proposed algorithm has been proven and demonstrated by means of numerical experiments. The algorithm has been also applied to a cosimulation of the power converter and the LHC main dipole circuit.
The proposed waveform relaxation algorithm will become an integral part of the STEAM framework and will allow for a more detailed analysis of protection circuits of superconducting magnets [1]. One application consists of cosimulation of a power converter controller, a circuit, and the a quenching magnet. Such a model would allow to study the influence of the quench initiation on the power converter’s response, and ensuing transients occurring in the circuit. A next step includes the analysis of the convergence of the electrical circuit with nonlinear elements and possibly partial differential equations.
Acknowledgment
The authors would like to thank Samer Yammine from CERN for a discussion on the LHC main dipole power converter structure and characteristics. This work has been partially supported by the Excellence Initiative of the German Federal and State Governments and the Graduate School of CE at TU Darmstadt.
References
 L. Bortot, M. Maciejewski, A. M. Prioli, Marco Fernandez Navarro, J. B. Ghini, B. Auchmann, and A. P. Verweij, “A consistent simulation of electrothermal transients in accelerator circuits,” IEEE Trans. Appl. Super., vol. 27, no. 4, Jun. 2016.
 I. Cortes Garcia, S. Schöps, L. Bortot, M. Maciejewski, M. Prioli, A. Fernandez Navarro, B. Auchmann, and A. Verweij, “Optimized field/circuit coupling for the simulation of quenches in superconducting magnets,” 2017.
 S. P. Panda, K. A. Salunkhe, and A. M. Kulkarni, “Experimental validation of waveform relaxation technique for power system controller testing,” Sadhana, vol. 40, no. 1, pp. 89–106, 2015.
 E. Lelarasmee, A. E. Ruehli, and A. L. SangiovanniVincentelli, “The waveform relaxation method for timedomain analysis of large scale integrated circuits,” IEEE Trans. Comput. Aided. Des. Integrated Circ. Syst., vol. 1, no. 3, pp. 131–145, 1982.
 J. K. White, F. Odeh, A. L. SangiovanniVincentelli, and A. E. Ruehli, “Waveform relaxation: Theory and practice,” Transactions of the Society for Computer Simulation, vol. 2, no. 1, pp. 95–133, 1985.
 K. Burrage, Parallel and sequential methods for ordinary differential equations. Oxford: Oxford University Press, 1995.
 C.W. Ho, A. E. Ruehli, and P. A. Brennan, “The modified nodal approach to network analysis,” IEEE Trans. Circ. Syst., vol. 22, no. 6, pp. 504–509, Jun. 1975.
 D. Estévez Schwarz and C. Tischendorf, “Mathematical problems in circuit simulation,” Math. Comput. Model. Dyn. Syst., vol. 7, no. 2, pp. 215–223, 2001.
 M. Günther and P. Rentrop, “Numerical simulation of electrical circuits,” GAMM, vol. 12, pp. 51–77, 2000.
 Z. Bartoszewski and M. Kwapisz, “On error estimates for waveform relaxation methods for delaydifferential equations,” SIAM J. Numer. Anal., vol. 38, no. 2, pp. 639–659, 2000.
 Z. Jackiewicz and M. Kwapisz, “Convergence of waveform relaxation methods for differentialalgebraic systems,” SIAM J. Numer. Anal., vol. 33, no. 6, pp. 2303–2317, Dec. 1996.
 A. Bartel, M. Brunk, and S. Schöps, “On the convergence rate of dynamic iteration for coupled problems with multiple subsystems,” J. Comput. Appl. Math., vol. 262, pp. 14–24, May 2014.
 E. Ravaioli, K. DahlerupPetersen, F. Formenti, J. Steckert, H. Thiesen, and A. Verweij, “Modeling of the voltage waves in the LHC main dipole circuits,” IEEE Trans. Appl. Super., vol. 22, no. 3, pp. 9 002 704–9 002 704, Jun. 2012.
 E. Ravaioli, K. DahlerupPetersen, F. Formenti, V. Montabonnet, M. Pojer, R. Schmidt, A. Siemko, M. S. Camillocci, J. Steckert, H. Thiesen, and A. Verweij, “Impact of the voltage transients after a fast power abort on the quench detection system in the LHC main dipole chain,” IEEE Trans. Appl. Super., vol. 22, no. 3, pp. 9 002 504–9 002 504, Jun. 2012.