Multistep Model Predictive Control for NPC Inverter Driving an Induction Machine
Overview
This RT Box model features a medium voltage three-level neutral-point-clamped voltage source inverter driving an induction machine. The squirrel cage induction machine has a rated power of 2 MVA. The objective of the model predictive current controller is to control the stator currents along their time-varying reference, by manipulating the switch position, while minimizing the switching effort. The predictive controller uses a prediction horizon of \(N_\mathrm{p}=5\). An efficient optimization algorithm is used to solve the underlying optimization problem.
The discretization step size and average execution times for this model are shown in Fig. 1. Note that the Controller execution time refers to the most loaded CPU core running the MPC algorithm with optimization enabled on an RT Box 2, 3 or 4.
Fig. 1 Performance overview for the execution on two RT Boxes
Requirements
To run this demo model, the following items are needed (available at www.plexim.com):
One PLECS RT Box 2, 3 or 4 to run the controller model and a second PLECS RT Box 1, 2, 3 or 4 to run the plant model
Follow the step-by-step instructions on configuring PLECS and the RT Box in the Quick Start guide of the RT Box Target Support Package Documentation.
Three 37 pin Sub-D cables to connect the boxes front-to-front.
Note
This demo model is designed for a two-RT Box setup, with an RT Box 1, 2, 3 or 4 running the Plant and an RT Box 2, 3 or 4 running the Controller. This configuration minimizes the execution time of each real-time target and allows a straightforward transition to HIL or RCP testing at a later stage. If only a single RT Box 2, 3, or 4 is available, refer to the corresponding PLECS model targeted for a one-RT Box setup. In this case, two 37-pin Sub-D cables are still required to connect the front-panel Analog Out to the Analog In interface, and the Digital Out to the Digital In interface.
Note
Model
The top-level schematic contains two subsystems representing the controller and the plant models, as shown in Fig. 2. Both subsystems are enabled for code generation from the Subsystem + Execution settings… context menu. This step is necessary to generate the model code for a subsystem via the PLECS Coder.
Fig. 2 Top-level schematic of the demo model
Plant
The schematic as seen in Fig. 3 shows a medium voltage (MV), three-level, three-phase, neutral-point-clamped (NPC) voltage source inverter driving a squirrel cage induction machine. The neutral point potential is fixed at zero and the dc-link voltage is assumed to be constant.
Fig. 3 Power circuit of the drive system
The three-level NPC is represented by three 3-Level Half Bridge (NPC) power modules. The pulse-width modulated (PWM) switching signals are obtained from the PWM Capture block of the PLECS RT Box library. Further details on the power module components and the sub-cycle averaging of PWM signals are described in [1]. The AC output current measurements are connected to Analog Out block from the PLECS RT Box library. The rotor angular position and rotational speed are converted by the Incremental Encoder block from the PLECS RT Box library, into digital orthogonal pulses, which can be measured outside of the subsystem.
Machine Parameters
The parameters of this model are borrowed from the case-study defined in [2]. The machine is a 3.3 kV, 50 Hz squirrel-cage induction machine rated at 2 MVA. The rated values of the machine and the drive parameters are provided in the tables below.
Parameter |
SI value |
|---|---|
Voltage |
3300 V |
Current |
356 A |
Real power |
1.587 MW |
Apparent power |
2.035 MVA |
Angular stator frequency |
\(2\pi 50\) rad/s |
Rotational speed |
596 rpm |
Air-gap torque |
26.2 kNm |
Parameter |
SI value |
pu value |
|---|---|---|
Stator resistance |
57.61 \(\mathrm{m}\Omega\) |
0.0108 |
Rotor resistance |
48.89 \(\mathrm{m}\Omega\) |
0.0091 |
Stator leakage inductance |
2.544 mH |
0.1493 |
Rotor leakage inductance |
1.881 mH |
0.1104 |
Main inductance |
40.01 mH |
2.349 |
Number of pole pairs |
5 |
|
dc-link voltage |
5.2 kV |
1.930 |
Controller
The controller schematic is shown in Fig. 4. The measurements of the AC currents are imported by Analog In block. The mechanical angular speed of the rotor is obtained from the Quadrature Encoder Counter block, which converts the orthogonal digital pulses.
Fig. 4 Predictive controller model of the drive system
This model uses model predictive control (MPC) scheme for the NPC inverter. The outer flux and speed control loops provide an input to the MPC controller.
The MPC controller predicts the stator currents at the next time step for all possible switch positions. For the prediction, the sampled stator currents are required, along with the rotor flux vector. The rotor flux vector cannot be measured directly. Therefore, this vector has to be estimated based on known quantities, i.e. the stator currents, stator voltages and rotor speed.
Control Problem
The outer flux and speed controllers provide the reference d- and q-component of the stator currents to the predictive current controller. The flux controller ensures that the machine is appropriately fluxed, while the speed controller ensures that the machine is close to the desired speed reference.
The predictive current controller can be split in two main parts:
the prediction of the future current trajectories based on a dynamic internal model of the drive system
the optimization problem to find the unique switching vector that minimizes the cost function solved with the modified sphere decoding algorithm.
Speed Controller
The speed value measured by the incremental encoder and the speed reference is provided as an input to the speed controller. The speed controller then regulates the electromagnetic torque setpoint, based on which the q-component of the stator current reference is calculated and provided to the predictive current controller.
Flux Controller
The flux controller ensures that the induction machine is appropriately fluxed and can achieve fast control of the flux magnitude. The error between the estimated rotor flux and its reference is provided as an input to the flux controller. The flux controller then manipulates the d-component of the stator currents.
Rotor Flux Estimation
Based on the measured stator currents, rotor speed and stator voltages, the rotor flux vector can be estimated. The stator voltage values can be reconstructed using the upper and lower dc-link voltages and the currently applied switch positions. In this demo the rotor flux estimation is just an implementation of the machine equations. Since an estimator has no corrective feedback, this rotor flux estimation scheme will not account for parameter variations, and parameter uncertainties.
Direct Model Predictive Current Control
MPC incorporates a dynamic model of the system to be controlled. This model is needed to predict a sequence of future system states and outputs based on the actual state and a sequence of manipulated variables. Based on a cost function the optimum future input sequence in order to drive the system towards a reference state is searched. The control objective is translated into a cost function, which maps the sequence of future states, outputs, and inputs into a scalar value. This scalar value is used to compare different possible future input values together and to chose the optimum one, i.e the solution with the lowest cost.
A commonly used MPC approach in power electronics is to directly manipulate the switch positions of the semiconductors in order to track a reference. In this demo the direct MPC formalism regulates the stator currents towards the reference values provided by the outer speed and flux controllers. In direct MPC the current controller and the modulator are grouped into one computational stage. Therefore, no modulator is needed in direct MPC.
The major disadvantage of direct MPC is the computational burden of solving the optimization problem especially for longer prediction horizons, since the number of possible switching sequences grows exponentially.
Tuning Variables
Contrary to classical control schemes, in direct MPC there are only two parameters that can be tuned:
the sampling interval \(T_\mathrm{s}\) and
the switching penalty \(\lambda_\mathrm{u}\)
The switching penalty is always a trade-off between tracking accuracy and switching effort for a given sampling interval. This weighting coefficient is used directly in the cost function as shown below. A lower bound on the sampling interval is mainly given by the computational power of the controller, and secondly the minimum on-time of the semiconductors.
Cost Function and Optimization Problem
The future values of the load currents are predicted for all the switching states generated by the inverter. For this purpose, it is necessary to measure the present load currents. After obtaining the predictions, a cost function \(g\) is evaluated for each switching state. The switching state that minimizes the cost function is selected and applied during the next sampling period.
The control requirements for the NPC inverter are:
Load current reference tracking
Reduction of the switching frequency to reduce switching losses.
These requirements can be formulated in the form of a cost function to be minimized. The cost function for the NPC inverter with a prediction horizon of 1 has the following composition:
where, the first two terms are the load current errors in orthogonal coordinates, where \(\mathrm{i_\alpha^p}\) and \(\mathrm{i_\beta^p}\) are the real and imaginary components of the predicted current vector \(\mathrm{i^p}\), respectively, and \(\mathrm{i_\alpha^*}\) and \(\mathrm{i_\beta^*}\) are the real and imaginary components of the reference current vector \(\mathrm{i^*}\).
The last term in Eq. (1) is the number of commutations \(\mathrm{n_c}\) required to switch from the present switching state to the switching state under evaluation. A switching state that implies fewer commutations of the power semiconductors will be preferred. In this manner, the use of this term will have a direct effect on the switching frequency of the converter. Weighting factor \(\mathrm{\lambda_n}\) handles the reduction of switching frequency within the cost function \(g\). A larger weighting factor value implies greater priority to that objective.
If the prediction horizon is increased, the cost function must be calculated for each step and finally summed up.
Computational Burden
The above optimization problem can be solved based on exhaustive search. In doing so, the cost function is evaluated for each possible switching combination. As the prediction horizon is enlarged and the number of decision variables is increased, the computational complexity grows exponentially. Therefore, exhaustive search is computationally feasible only for very short predictions horizons. However, a far more efficient method based on branch-and-bound for finding the optimal switching sequence can be used. The second method is an adapted version of a so-called sphere decoding algorithm.
Modified Sphere Decoding Algorithm
The modified sphere decoding algorithm is based on the branch-and-bound approach. It is dentified to be a very efficient algorithm to find the optimal switching sequence \(U_{\mathrm{opt}}\) for direct model predictive control optimization problems with long prediction horizons, by only exploring a very small subset of the search tree. Although only a small part of the possible nodes is calculated, the algorithm still returns the optimal solution in each sampling period. For completeness, the pseudo-code of the modified sphere decoding algorithm introduced in [2] is given in Fig. 5.
Fig. 5 Modified sphere decoding algorithm
To invoke the search algorithm, it must be initialized. For this, first the unconstrained optimal solution \(\bar{U}_{\mathrm{unc}}\) can be calculated explicitly, as shown in [2]. This unconstrained solution neglects the fact that a valid switching vector can only contain discrete values as given in \(U_{\mathrm{set}}\) [1]. Direct componentwise rounding of the unconstrained solution towards the nearest integer value given in \(U_{\mathrm{set}}\)
\[U_{\mathrm{ini}} = round(U_{\mathrm{unc}})\]
leads to a possible sub-optimal candidate solution. However, this sub-optimal solution can be used as an initial guess to invoke the optimization routine. A sphere centered around \(\bar{U}_{\mathrm{unc}}\) is spanned and the radius \(\rho_{\mathrm{ini}}\) is set in such a way that the initial solution \(U_{\mathrm{ini}}\) just lies within the sphere. This ensures that at least one valid solution will be found. As shown in the pseudo-code of the modified sphere decoding algorithm, a possible switching vector \(U\) is built up componentwise in a recursive way as long as it lies within a sphere described by the radius \(\rho\). This is the case, if the squared distance between the unconstrained solution vector and the current switching vector candidate \(U_{1:i}\) is less than the current squared distance \(\rho^2\). If the candidate vector is not inside the current sphere, the search in this branch of the tree is aborted immediately since it will lead to a sub-optimal solution. Once the vector has full length, e.g. \(i = 3N_\mathrm{p}\), a better solution was found and the radius of the sphere is reduced to the new squared distance between the unconstrained solution and the new best solution. By doing so, the sphere radius is tightened every time a better switching vector is found, and at the end of the routine only the single switching vector minimizing the cost-function lies within the sphere. A complete mathematical derivation of the control problem is given in [2].
The sphere decoding algorithm is invoked as
and returns the optimal switching vector \(U_{\mathrm{opt}}\) that is applied to the gate drivers.
Simulation
This model can run both in offline mode on a computer, and in real-time mode on the PLECS RT Box. For the real-time operation, two RT Boxes (referred to as “Plant” and “Controller”) need to be set up as demonstrated in Fig. 6 using three DB37 cables.
Please follow the instructions below to run a real-time simulation on two RT Boxes (one RT Box 1, 2, 3 or 4 for the “Plant” and one RT Box 2, 3 or 4 for the “Controller”):
Connect:
Analog Out interface of the “Plant” RT Box to the Analog In interface of the “Controller” RT Box
Digital Out interface of the “Plant” RT Box to the Digital In interface of the “Controller” RT Box
Digital In interface of the “Plant” RT Box to the Digital Out interface of the “Controller” RT Box.
From the System tab of the Coder options… window, select the “Plant” and Build it onto the “Plant” RT Box. Then, select “Controller” and Build it onto the “Controller” RT Box.
Once both models are uploaded, from the External Mode tab of the Coder options… window, Connect to both RT Boxes and Activate autotriggering.
Fig. 6 Hardware configuration for the real-time operation of the demo model
Key waveforms can be observed in the different Scopes individually placed in the “Plant” and “Controller” subsystems. The optimization of the predictive current controller can be enabled by toggling the “Enable Optimization” manual switch inside the controller subsystem from 0 to 1. If the optimization is disabled, the sphere decoding algorithm is not invoked and a sub-optimal explicit switching vector is applied to the inverter.
The Scope “Explored nodes” shows the number of switching combinations the optimization algorithm has explored before the optimal solution was found. This Scope is only relevant once the optimization has been enabled as explaind above.
The Manual switch to enable a load torque step in the plant subsystem is inlined and can therefore be toggled during a real-time simulation. When toggling the manual switch a load step of 0.5 p.u is applied to the drive system.
The main simulation results are shown in Fig. 7.
Fig. 7 Steady-state voltage and current waveforms for a prediction horizon of \(N_\mathrm{p} = 5\) and a load torque of 0.5 p.u. The cost function weighting factor was chosen to achieve a switching frequency around 300 Hz
In addition, the prediction horizon can be changed prior to an offline and real-time simulation. From the Model initialization commands window of Simulation Parameters… + Initializations tab from the Simulation menu, change the value of Par.Ctr.Np to a value between 1 and 5, to choose the desired prediction horizon. You must also configure the prediction horizon in the “MPC Controller” C-Script in the controller subsystem accordingly. Open the C-Script dialog and navigate to the Code declarations section in the Code tab. Change the value of Np to the same value as in the Model initialization commands. You can see the effect of different prediction horizons in the “Explored nodes” scope. The lower the prediction horizon is, the less nodes have to be explored by the sphere decoding algorithm. Please note, that the model must be reloaded on both RT boxes after changing the prediction horizon:
Open the web interface of each RT Box by clicking on the
icon under the Target or the External Mode tabs of the Coder Options dialog.Stop the execution of the real-time simulation by clicking on the Stop button at the bottom of the Application tab in the web interface.
Re-build the “Plant” and “Controller” subsystems on the respective RT Boxes as described in Section Simulation.
Conclusion
This model demonstrates a medium voltage three-level neutral-point-clamped voltage source inverter driving an induction machine. The drive system is closed-loop controlled by means of outer flux and speed controllers and an inner predictive current controller. The model can run both in offline and real-time operation for Hardware-in-the-loop testing and rapid control prototyping.
Bibliography
[1]
J. Allmeling, and N. Felderer, “Sub cycle average models with integrated diodes for real-time simulation of power converters,” IEEE Southern Power Electronics Conference (SPEC), 2017.
[2]
T. Geyer and D. E. Quevedo, “Multistep Finite Control Set Model Predictive Control for Power Electronics,” in IEEE Transactions on Power Electronics, vol. 29, no. 12, pp. 6836-6846, Dec. 2014, doi: 10.1109/TPEL.2014.2306939.
[3]
L. Capponi, “Direct Model Predictive Control Using Model Based Design,” IEEE International Conference on Predictive Control of Electrical Drives and Power Electronics (PRECEDE) , 2021.