First application of data assimilation-based control to fusion plasma

Magnetic fusion plasmas, which are complex systems comprising numerous interacting elements, have large uncertainties. Therefore, future fusion reactors require prediction-based advanced control systems with an adaptive system model and control estimation robust to uncertainties in the model and observations. To address this challenge, we introduced a control approach based on data assimilation (DA), which describes the system model adaptation and control estimation based on the state probability distribution. The first implementation of a DA-based control system was achieved at the Large Helical Device to control the high temperature plasma. The experimental results indicate that the control system enhanced the predictive capability using real-time observations and adjusted the electron cyclotron heating power for a target temperature. The DA-based control system provides a flexible platform for advanced control in future fusion reactors.


Methods
Figure 1 shows a conceptual diagram of the DA-based control system.This system has a simple structure comprising a DA system, heating and fuel supply equipment, and measurement devices.This section describes the DA-based control system constructed for the LHD plasma control.

DA framework for adaptive predictive control
The DA system, which is the core component of the control system, computes the probability distribution of the state vector (state distribution), which is defined at time t as x t = (x T t , u T t ) T , where vector xt is the part of x that includes the system state and model parameters, and vector u is the control input that determines the time evolu- tion of the system (see Table 1 for an example).The superscript "T" denotes the transpose of a vector and t z is the time interval for estimating the control inputs.The DA system optimizes x based on the observed information y (adaptation) and estimates u based on the predicted state distribution and target state z (predictive control).Consider a control problem in which the control input is adjusted at each time interval t z and the system state is observed as the vector y at each time interval t y .For simplicity, assume that t y = n t z ( n ∈ N ) and use the time notation (i, j) = t i,j = t 0 + i�t y + j�t z and t i = t i,0 , where t 0 denotes the initial time.The DA framework for adaptive predictive control 18 is based on the following state-space model: Equation (1) represents the system model that describes the time evolution of the system x (i,j) → x (i,j+1) , consid- ering the effect of system noise v (i,j+1) .Assume that the value of u t is constant in the prediction interval t z , i.e., Here, the system noise for the control input v u (i,j+1) is added to u (i,j) before calculating the time-evolution.Equa- tions ( 2)-( 4) represent the relationship between the state vector x (i,j) and vectors z (i,j) , u * (i,j) , and y i based on the noises w z (i,j) , w u (i,j) , and w y i , respectively.The matrices H z , H u , and H y are linear operators for projecting the state vector in each corresponding space.The system noise v (i,j+1) is assumed to follow a Gaussian distribution with zero mean and covariance matrix Q (i,j+1) , i.e., v (i,j+1) ∼ N(0, Q (i,j+1) ) .Similarly, w z (i,j) , w u (i,j) , and w y i are assumed to follow the probability distributions N(0, R z (i,j) ) , N(0, R u (i,j) ) , and N(0, R y i ) , respectively.The priority of the controlled and observed variables can be adjusted using w z (i,j) and w y i .In the DA system, the state distribution is approximated by an ensemble using parallel computing.Each ensemble member represents a simulation with slightly different conditions (e.g., initial conditions, model parameters, and control inputs).In this experiment, considering the time required for computation and communication, the prediction and control estimations were performed up to t y ahead of the real time at n = 3 ( t y = 3 t z ).The following are the computational steps for the adaptive predictive control: • y-filter here the subscript t 1 : t 2 denotes all the timings in [t 1 , t 2 ] .Given the distribution p(x (i,j) | y 0:i−1 , u * (0,1):(i,j) ) , the prediction step computes p(x (i,j+1) | y 0:i−1 , u * (0,1):(i,j) ) , the distribution t z ahead, based on the system model.This step can be performed by computing the time evolution of each ensemble member approximating the state distribution.
The z-and u-filter are the control estimation steps, and the y-filter is the adaptation step.These three steps are based on the Bayesian filter, which calculates the conditional probability distribution p(x | ξ) using the distribu- tion p(x) and the model representing the relationship between x and the imposed information ξ , ξ = h(x) + w .Here, h represents the relationship between ξ and x and w represents the associated noise.Generally, ξ denotes the observation data, and the Bayesian filter can assimilate the observed information into the state distribution.Kalman filter (EnKF) 19 and particle filter (PF) 20 are typical computational methods for the Bayesian filter that use an ensemble.
The z-filter step estimates the control input by assimilating the target z (i,j+1) into the predicted distribution based on the Bayesian filter using Eq. ( 2), and performing marginalization (1) We used EnKF as the computational method for the Bayesian filter.This marginalization can be performed by removing x from the ensemble of the distribution p(x (i,j+1) | y 0:i−1 , u * (0,1):(i,j) , z (i,j+1) ) .The control input u * (i,j+1) is obtained as the expected value of the z-filtered distribution.The u-filter is applied by assimilating the estimated control input u * (i,j+1) into the predicted distribution using the Bayesian filter based on Eq. ( 3).The control estimation process can proceed by repeating the following three steps: prediction→z-filter→u-filter.
The y-filter corresponds to the adaptive process and suppresses the uncertainties in the system model by assimilating the observations into the state distribution.Because the observation times t i differ from the latest u-filtered distribution t i+1,0 , the y-filter is executed by assimilating y i into the joint distribution of the state vectors at two time points.The ensemble approximating the joint distribution p(x (i+1,0) , x (i,0) | y 0:i−1 , u * (0,1):(i+1,0) ) can be obtained by concatenating the u-filtered ensemble at t i+1,0 with the stored ensemble at t i .The filtered distribu- tion p(x (i+1,0) , x (i,0) | y 0:i , u * (0,1):(i+1,0) ) can be calculated by assimilating y i with the concatenated ensemble using Eq. ( 4).The ensemble of p(x (i+1,0) | y 0:i , u * (0,1):(i+1,0) ) is obtained by marginalizing x (i,0) .This DA framework can be used to construct an adaptive predictive control algorithm without overlapping the prediction intervals.The control procedure for this experiment is summarized in Fig. 2. In t i < t < t i+1 , while the real system evolves in time with the inputs u * (i,1) , u * (i,2) , and u * (i+1,0) , the predictions and control estimates are performed from t i+1 to t i+2 in the DA system, as shown in Fig. 2a.At t = t i+1 , the system state is observed as y i+1 , which is assimilated into the latest u-filtered distribution (adaptation), as shown in Fig. 2b.The loops of these processes illustrated in Fig. 2 provide an adaptive predictive control.
The proposed control system for LHD comprises a DA system "ASTI" (Assimilation System for Toroidal plasma Integrated simulation) 21 , which uses an integrated transport simulation code TASK3D 22,23 as the system model for helical fusion plasmas and implements the EnKF and PF as DA techniques.TASK3D computes the heat and particle transport in a toroidal fusion plasma as a one-dimensional (1D) problem for the normalized minor radius.In this study, ASTI computed 256 ensemble members (TASK3D simulations) on a vector computer (NEC SX-Aurora TSUBASA, 16VE) connected to the LHD experimental system, and EnKF was used to perform the Bayesian filters.The computer had a maximum of 128 parallel processes, each of which was responsible for two ensemble members, thereby computing the 256 ensemble members.The control performance reached its maximum with 200 ensemble members in numerical experiments 18 to control the virtual LHD plasma generated by TASK3D using noise values similar to those in the LHD experiments.Although the behavior of the virtual plasma is different from that of the actual plasma, the system model and state vector are almost the same as in the experiment.Therefore, the effect of the number of ensemble members on the control performance is considered to be sufficiently small in the LHD experiment.

System model for LHD plasmas
We employed the TASK3D code 22,24 as the time evolution model f in Eq. (1).TASK3D is an integrated simulation code for helical fusion plasmas that solves the 1D diffusive transport problem in the radial ρ-direction.The parameter ρ is given by the magnetic flux surface, where 0 and 1 correspond to the center and edge of the plasma, (12)    www.nature.com/scientificreports/respectively.We solved the heat transport equations for each electron and ion species in the LHD experiment.
Assume that the electron and ion density profiles are similar, i.e., n = n e = n i .The radial profile in TASK3D was defined at 60 grid points.The geometric parameters required in the 1D transport simulation were evaluated by an equilibrium magnetic field precomputed using the VMEC code 25 for the typical LHD magnetic configuration.The major radius of the magnetic axis in vacuum was 3.6 m and the magnetic-field strength at the plasma center was 2.75 T. In TASK3D, changes in the geometric parameters were not considered in the time evolution calculation due to the computational cost.
For the electron and ion thermal diffusivities, the following gyro-Bohm models were employed: where B, ρ i , and a are the magnetic field strength, ion Lar- mor radius, and plasma minor radius, respectively.We set C gB e = 1.5 and C gB i = 0.1 as reasonable values based on previous studies 24,26 .The heating power source comprised the externally applied ECH, power exchange between species, and loss term by interaction with neutrals.ECH only contributes to the heating term of the electron and the following ECH model was employed in the experiments: where µ ECH = 0.15 and σ ECH = 0.04 are obtained from the ray-tracing calculations 27 .The coefficient A is deter- mined from the total ECH input power, which is given by where V is the plasma volume and V ′ = dV /dρ.

ECH control experiment
To demonstrate real-time adaptation and control estimation using the developed system, we considered a simple LHD experiment to control the central electron temperature by adjusting the ECH power 28,29 .Table 1 lists the state ( x and u ), target state ( z ), and observation variables ( y ) used in the experiment.The radial profiles of the state variables were defined on 11 grid points ( ρ = 0, 0.1, 0.2, . . ., 1 ) in the state vector, and B-spline interpolation 30 was used to transform the radial profiles to those on the TASK3D grid (60 grid points) 21 .The factors for the thermal diffusivities c e and c i were introduced to optimize χ e and χ i with large uncertainties.In the prediction step, c e χ e and c i χ i were used instead of χ e and χ i in the TASK3D simulation.The ECH used the four gyrotrons: 330, 400, 600, and 600 kW for 5.5 s; thus, ASTI controlled the injection power by selecting a subset of these gyrotrons by switching the individual gyrotrons on or off.ASTI sends the estimated control signals to the gyrotrons every t z = 0.1 s and receives the observed radial profiles of the electron temperature and density every t y = 0.3 s from the real-time Thomson scattering measurement system.
To ensure that the behavior of the system model resembles that of the real system at the beginning of control, observation assimilation was performed with a fixed heating power of 731 kW during the first phase t < 2.1 s.ASTI assimilates the observations up to 1.8 s, after which the control estimation begins.Subsequently, the ECH power control commences at 2.1 s to produce the target electron temperature (4 keV) at 3.9 s and maintain it.The electron density at the plasma center was maintained at 1.5 × 10 19 m −3 using the PID control to focus on the temperature control.
The covariance matrices for the noises Q (i,j+1) , R z (i,j) , R u (i,j) , and R y i are the key parameters affecting the overall control performance.In addition to them, the ensemble mean x(0,0) and covariance matrix V (0,0) are required to generate the initial ensemble members.Diagonal matrices were used for these covariance matrices, i.e., the covariance component of the matrices was not considered.
The covariance matrix for the system noise, Q (i,j+1) , controls the uncertainty of the system state, including the model parameters and control inputs.System noise was added before each prediction step, and x was assigned a slightly larger value of noise after processing by the y-filter to prevent the distribution from shrinking.The standard deviation of the system noise was fixed at the values listed in Table 2.The values of the standard deviation for T e , T i , c e , and c i were determined based on a previous study on data assimilation for the LHD plasmas 18,21 .The noises for n e were set to large values because ASTI did not compute the density transport.In control estimation, system noise for u is an important parameter determining the range of control inputs considered in a single control estimation and the change rate of u .In this experiment, the standard deviation of the noise for P ECH was set to a sufficiently large value to track the rate of change of the target temperature.
The covariance matrix R z (i,j) affects the performance of the z-filter and determines the proximity of the system state to the target state after the z-filter step.The diagonal components of R z (i,j) were determined at each z-filtering step as follows: where V (i,j) denotes the covariance matrix of the ensemble that approximates the predicted distribution at t i,j and r z = 0.5 18 .The subscript ( ) ll and superscript T denote the l-th diagonal component and matrix transposition, respectively.
The covariance matrix R u (i,j) considers the uncertainty in the control input.In the experiment, the standard deviations of the control input noise were set to a sufficiently small value of 0.05 MW for P ECH .
The covariance matrix R y i affects the performance of the y-filter and determines the effect of the observed information on the state distribution.The standard deviation of the observation noise is assumed to be proportional to the difference between the observation and mean of the state distribution 18 , where r y = 0.8 , x(i,0) is the mean of the ensemble approximating p(x (i,0) | y 0:i−1 , u * (0,1):(i,0) ) , and ( ) l represents the l-th element of the vector.This assumption prevents the control instability caused by overfitting of the system model to noisy observations and maintains the variance of the y-filtered ensemble at an adequate magnitude.
The initial ensemble means of T e and T i are set to the steady-state radial profiles calculated by the TASK3D simulation for the initial ECH.The initial mean profile of n is set to n(ρ) = 1.0 − 0.8ρ 8 [ ×10 19 m −3 ] and those of c e and c i are set to 1.The initial ensemble variances are assigned the values listed in Table 2.

Real-time observation and control system
The DA system ASTI was connected to the real-time measurement system and ECH system 28,29 via the real-time communication system in LHD 31 .The temperature and density profiles along with their measurement errors obtained from the Thomson scattering measurement system 32,33 were sent from the Thomson data analysis PC to the vector engine server (SX-Aurora TSUBASA, NEC Inc.) via socket communication.The digitizers of the highrepetition-rate Thomson scattering system 34 were used for real-time measurements (10 Hz).In the vector engine server, ASTI considers the time delay relative to real time, which is tens of milliseconds.This is attributed to the temperature computation, density analysis, and communication delay.The electron density and temperature were observed at 144 radial points, of which 66 data points were available for real-time measurement.After removing the obvious outliers based on measurement errors, a mapping model was used to transform the observed data in the real space coordinate (major radius R) into a form suitable for DA (profiles on 11 points of the ρ-coordinate).
The mapping model acted as coordinate transformation, outlier removal, data smoothing, and data point extraction for DA, providing the observation vector y to ASTI.The mapping model was based on a multilayer perceptron with five hidden layers containing 500 units each using ReLu as the activation function.The mapping model was trained using data on the LHD experimental database (profiles observed by the Thomson scattering system and the mapped data obtained from the equilibrium calculations).The training data comprised 58 discharges (7769 time points) of the ECH plasma, of which a quarter were used as the test data.For all discharges in the training data, equilibrium magnetic fields were calculated and used to evaluate the density and temperature profiles in the ρ-coordinate for the teacher data.Before training, the mapped data were fitted with even-order polynomials up to the 8th order to smooth the profiles and extract the data points required for DA.An even function was used for fitting because the Neumann boundary condition ∂/∂ρ = 0 was imposed on ρ = 0 in TASK3D.The variations of the magnetic surfaces were considered in this mapping model, and the observations in the ρ-coordinate were provided to ASTI.Since TASK3D did not solve for the magnetic equilibrium, errors occurred in the geometric parameters required for the 1D transport simulation (e.g., dρ/dR ).In this experiment, Table 2. Standard deviations of the initial state distribution ( σ init ) and system noise ( σ Q ).The values in parentheses represent the additional system noise added after the y-filter.The values with% as the unit represent the rate for determining the standard deviation in proportion to the state distribution mean.www.nature.com/scientificreports/these errors were considered small for helical devices and allowed.However, the real-time computation of the magnetic field, which is essential for tokamak control, will be addressed in future studies.
After the prediction and control estimation in ASTI were completed, the generated ECH control signals were sent to the computer (Jetson AGX Orin, NVIDIA Inc.) via socket communication.The digital signal for ECH control, which was converted to an analog signal using a D/A converter (EVAL-AD5686RSDZ, Analog Device Inc.) connected to a computer, was input to the ECH system to control the ON/OFF state of the ECH injection.

Results
Figure 3 shows the control results of the LHD experiment and Fig. 4 shows the prediction error (the absolute difference between "Prediction" and "Observation" in Fig. 3).The electron temperature increases to that of the target (4 keV) and is maintained beyond 3.9 s.The prediction error, shown in Fig. 4, increases to 2 keV in the transient section (2.1 s ≤ t < 3.9 s) and decreases to approximately 0.2 keV in the steady-state section (3.9 s ≤ t ).These results demonstrate the effectiveness of real-time adaptation and control estimation using the DA-based  Prediction error defined as the absolute difference between "Prediction" and "Observation" in Fig. 3. control system.The prediction error in the transient section occurred mainly because the gyro-Bohm model overestimated the thermal diffusivities as the electron temperature increased.At the beginning of the control (23 s), the c e optimization lagged behind the variation in the actual thermal diffusivities.However, before 4 s, the gap between the model and the actual system was bridged, which improved the prediction performance.The prediction error also depends on the Thomson measurement error of 0.2-0.4keV and accuracy of the mapping model.
In addition, the discrete nature of the control input and expansion of the u-filtered distribution from 3.4 s onward also affect the control accuracy.This expansion is attributed to the increase in ensemble members with small c e (< 0.3) near the center as shown in Fig. 3c, which increases the range of electron temperature variation.Because the state distribution is approximated by finite ensemble members, a large uncertainty in the predicted distribution makes it difficult to identify the relationship between the control inputs and the system state and may reduce the control accuracy.The expansion of the u-filtered distribution at 3.3 s was triggered by adding system noise with the fixed standard deviation of 0.1 to c e , despite the small mean of c e near the center.This expansion can be avoided by introducing a method to dynamically adjust the system noise by considering the characteristics of the system model.The probability distribution of model parameters should be adjusted by considering the physical characteristics and control accuracy.This issue will be addressed in a subsequent control experiment.

Change of state distribution in the control process
The ensemble members of the state distributions at 3.0 s are shown in Fig. 5. Figure 5a shows the z-filter process in which the target temperature is assimilated into the predicted ensemble.The z-filtered ensemble approximates the distribution p(x 3.0 | y 0.3:2.4 , u * 2.1:2.9 , z 3.0 ) , where the subscript number indicates the time (s) and t 1 : t 2 denotes the values in [t 1 , t 2 ] .As shown in Fig. 5b, the u-filtered ensemble is computed after the control input u * 3.0 is obtained as the expected value of this distribution.The distribution approximated by the u-filtered ensemble corresponds to the predicted distribution with zero or small uncertainty in the control input.The remaining uncertainties originate from the uncertainties in the model parameters and system state before time evolution calculation.The u-filtered ensemble is modified by assimilating the observation y 2.7 (y-filter), as shown in Fig. 5c.This filter optimizes the state variables, including the model parameters, to adapt the system model to the real system.The control process proceeds to the subsequent prediction from the y-filtered ensemble.
The filtering steps were performed using EnKF in the simple control experiment.Strongly nonlinear relationships among the state variables and a significantly non-Gaussian state distribution (e.g., a multimodal distribution) necessitate the use of other filters (e.g., PF) and another definition of u * (e.g., mode).

Real-time adaptation to the real system
The upper panels of Fig. 6 show the radial profiles of the predicted electron temperature (u-filtered distribution) and observed profiles at three different points in time.The radial profiles of the electron thermal diffusivity are shown in the lower panels.The thermal diffusivity can be estimated from the observed density and temperature profiles as where represents the magnetic flux surface average.This estimate is less valid near the center where the temperature gradient is smaller.
After the first phase (observation assimilation), the temperature profile was predicted with high accuracy at the beginning of control estimation (t = 1.8 s), as shown in Fig. 6a,d.In the transient phase (2.1 s ≤ t < 3.9 s), insufficient modeling of the electron thermal diffusivity (Fig. 6e) underestimated the profile around the plasma center, as shown in Fig. 6b.These prediction errors contribute to the control error in the transient phase, as shown in Fig. 3.However, in the steady-state phase (3.9 s ≤ t ), the temperature profile was predicted with high accuracy by optimizing the thermal diffusivity model, as shown in Fig. 6c,f.
These results demonstrate the real-time adaptation of the integrated transport simulation code to the actual LHD plasma.In the transient phase, changes in the circumstances and the discrepancy between the system model and real system degrade the prediction performance.Furthermore, in the control procedure used in this experiment, the observations obtained at time t i were reflected in the control estimation after t i + t y , which caused a delay in adaptation, especially at the beginning of the control (2.1 s−).This delay can be a critical issue for some control problems.Therefore, to achieve a stable adaptive predictive control system, it is important to design the control procedure, t y , and the observation noise to minimize this delay.
Reduction of the control error in the transient phase makes it necessary to improve the assumed transport model, adjust the adaptive capacity by controlling the uncertainty in model parameters, and minimize the adaptation delay.These limitations will be addressed in the future studies aimed at building a stable control system.The DA framework can be applied to construct an advanced control system that can optimize the timing, physical quantities, and measurement positions.For instance, it would be possible to develop a system that allows measurement devices to observe the system state depending on the uncertainty of the model parameters and system state.

Discussion
A control system based on DA was developed to control LHD plasmas.A simple experiment to control the electron temperature demonstrated adaptive predictive control using a nonlinear system model (integrated transport model) with large uncertainties.The proposed control framework simplifies the structure of the control system to accommodate complex control problems characterized by a large number of variables.
The adaptive capacity of a control system is limited by the assumed models and frequency of observations 18 .Therefore, the sophistication of the model, computational cost, and amount of information contained in the observations must be adjusted depending on the target control problem.Similar to other control systems, observation errors significantly affect the performance of the proposed DA-based control system.However, this system can explicitly account for the uncertainties in observations based on the observation noise.Robust and stable control can be achieved by exploiting this uncertainty; however, methods to adjust the uncertainties in observations and model parameters according to the situation should be developed.
Advanced control systems constructed for future fusion reactors would require models based on data-driven approaches, such as surrogate models [35][36][37] , stability indicators 15,38,39 , and tomography methods 40,41 .To this end, the DA-based control system, which expresses the system state as a probability distribution, is highly compatible with other data-driven approaches and provides a flexible control platform that links the physics-based and data-driven models.
The results of the first application of DA-based control system in LHD do not exceed the control performance of conventional controllers (e.g., PID controller).However, the proposed approach has high potential in the control of future fusion reactors under large uncertainties because it can estimate and control the internal state from limited observed information and naturally handle various observed and control variables (e.g., temperature and density profiles) and plasma responses with different time scales in a single framework.Thus, the success of this demonstration paves the way for addressing challenging control problems in fusion plasma, such as multivariate control, radial profile control, electric field bifurcation control, and avoidance of terminating events using relevant alarm rates.DA-based control may become a key technology for solving challenging control problems in future fusion reactors and other complex systems by augmenting numerical models with real systems.Beginning in 2024, our roadmap includes conducting more sophisticated control experiments with this control system by incorporating additional measurement, heating, and fuel-supply devices into the system.Additionally, we are actively working on extending our control capabilities to include Tokamak control.

Figure 1 .
Figure 1.Overview of the DA-based control system for fusion plasmas.(a) The state distribution in the DA system (ASTI) is expressed as an ensemble.The parallel time evolution of each ensemble member approximates that of the state distribution.(b) The control estimation is performed by assimilating the target state into the predicted state distribution.

Figure 2 .
Figure 2. Control procedure of the LHD experiment ( t y = 3 t z ).The control is performed by repeating the control estimation ( a ) and DA of the observation into the latest u-filtered distribution ( b).

Figure 3 .
Figure 3. Results of a control experiment (#186500).( a ) Control result of T e at the plasma center.The plotted values labeled "Prediction" correspond to the expected values of the predicted distribution for t ≤ 2.1 s and those of the u-filtered distributions for t > 2.1 s.The shaded areas represent a single standard deviation of the distributions.The plotted values labeled "Observation" are those obtained by the mapping model from the Thomson scattering measurements.( b ) ECH power adjusted using ASTI.( c ) and ( d ) Distribution (expected value and one standard of ( c e ) at ρ = 0.2 and 0.4 used in the prediction step ("Prediction") and y-filtered distribution ("y-filtered").

Figure 4 .
Figure 4. Prediction error defined as the absolute difference between "Prediction" and "Observation" in Fig.3.

Figure 5 .
Figure 5. Variations in the ensembles that approximate the state distributions at t = 3.0 s in the control experiment.(a) Estimation of the control input (z-filter).(b) Predicted distribution modified by the estimated control input (u-filter).(c) Adaptation to the actual LHD plasma by assimilating the observation at t = 2.7 s (y-filter).

Figure 6 .
Figure 6.Time evolution of the radial profiles of electron temperature and thermal diffusivity.(a-c) Radial profiles of the predicted electron temperature corresponding to Fig. 3a and the observations (mapped to the ρ-coordinate) at three points in time.(d-f) Radial profiles of electron thermal diffusivity calculated in the prediction step ("Prediction") and the values estimated using Eq.(17) from the observed temperature profile at the corresponding time points ("Estimated value").The shaded areas around the radial profiles represent a single standard deviation. https://doi.org/10.1038/s41598-023-49432-3

Table 1 .
State, target, and observation variables with their dimensions in the vectors ( M i ).