esa-sceva/satcom-chunk-collection huggingface dataset:
# Prediction of Corroded Pipeline Failure Pressure Based on Empirical Knowledge and Machine Learning
###### Abstract  
This paper presents a novel approach for predicting the failure pressure of corroded pipelines by integrating empirical formulas into the loss function of a neural network-based prediction model. Traditional empirical formulas, such as ASME-B31G, DNV RP-F101, and PCORRC, have been widely used for their simplicity but often suffer from significant prediction errors due to the complex interactions between defect parameters and material properties. In contrast, artificial neural networks (ANNs) offer more accurate predictions but require substantial training data. To address these limitations, we propose an integrated loss function that combines the strengths of empirical formulas and the powerful fitting capabilities of ANNs. The proposed loss function incorporates an additional defect factor term predicted by the neural network to compensate for errors caused by varying defect conditions, thereby enhancing the model's adaptability and accuracy. The model is trained using a diverse dataset of 60 burst test results from various literature sources, covering a wide range of corrosion scenarios. The results demonstrate that the proposed method significantly improves prediction accuracy compared to traditional empirical formulas and ANN models trained with standard loss functions. The proposed approach achieves a mean absolute percentage error (MAPE) of 2.52%, a root mean square error (RMSE) of 0.39 MPa, and a coefficient of determination (R\({}^{2}\)) of 0.9886 on the validation set. This study highlights the effectiveness of integrating empirical knowledge with data-driven models and provides a robust and accurate solution for predicting the failure pressure of corroded pipelines, contributing to enhanced pipeline integrity assessment and safety management.  
corroded pipeline; failure pressure prediction; empirical formula; loss function +
Footnote †: journal: Applied Sciences
# On Efficiency of Two-Degree-of-Freedom Galloping Energy Harvesters with Two Transducers
###### Abstract  
This paper examines the energy efficiency of three variations of the two-degree-of-freedom transverse galloping energy harvester. These variants differ in the number and placement of electromechanical transducers. By utilizing the harmonic balance method, the limit cycles of mathematical models of the devices were determined. Analytical expressions derived from the models were then used to formulate the efficiency of the systems. It was demonstrated that efficiency depends on flow speed and can be comprehensively characterized by the following criteria parameters: peak efficiency, denoting the maximum efficiency of the system, and high-efficiency bandwidth, which describes the range of flow speeds within which the efficiency remains at no less than 90% of peak efficiency. The values of these parameters are heavily reliant on two other parameters: the speed at which the system achieves peak efficiency, referred to as the nominal speed, and also the flow speed at which the system undergoes Hopf bifurcation, namely the critical speed. Comparative analysis revealed that only the device equipped with two electromechanical transducers can potentially outperform a simple one-degree-of-freedom system. For selected parameters, this gain reached nearly 10%.  
energy harvesting; transverse galloping; harmonic balance method; nonlinear dynamics; flow-induced vibration +
Footnote †: journal: Applied Mechanics Letters  
1
Footnote 1: email: {[[https://doi.org/10.3390/app14135427](https://doi.org/10.3390/app14135427)]([https://doi.org/10.3390/app14135427](https://doi.org/10.3390/app14135427))}
## 3 Proposed Design
### Numerical Analysis
It can be concluded that if there is a need to send 6 frames, an A-MPDU can sent it in one channel access with a minimum delay as in the cost of one frame it is sending five more frames. In contrast, the MPDU takes at least six channel attempts and the delay is significantly high. We also incorporate this analysis in our simulations. Figure 9 presents the MT analysis for the same scenario mentioned in Figure 8. It is shown that MT increases with the increase of aggregation level like from two frames in an A-MPDU to maximum of six frames. Therefore, it is concluded that, frame aggregation mechanism enhances the packet delivery which is helpful to reduce the delay and provision of time bounded services in WBASNs. The maximum throughput (MT) is considered as the ratio of payload size (bytes) to the transmission delay for specific payload size and is computed as:  
\[\text{MT}=\frac{\text{(x)}\times 8}{\text{transmission delay}} \tag{13}\]  
Figure 8: Delay comparison between MPDU and A-MPDU. Figure 9: Maximum throughput analysis.
## 3 Construction of LPTN-Based VML
### Calculation and Classification of TRs
#### 3.3.1 Constant TRs  
The constant TRs in the bearing LPTN model primarily include the radial conduction TRs of the inner/outer rings, conduction TRs of rolling elements, conduction TRs of the rotating shaft and bearing chamber, and contact TRs at solid-solid interfaces. Due to the negligible influence of external loads on contact TRs, these terms are omitted from computational considerations. Additionally, the natural convection heat transfer coefficient between the bearing chamber's outer surface and ambient air is treated as a constant. Mathematical formulations for these constant TRs are summarized in Table 2[49].
## 2 Problem Definition and Proposed Models  
The considered problem is to check the feasibility of machine learning for accurately estimating the CPU consumption of different virtual network functions deployed in SDN  
Figure 2: Basic operation of SDN-with-knowledge.  
environment, when trained with only traffic features as an input. The CPU needs of a particular VNF depending upon many known and unknown interacting parameters, which is the key motivation factor of using ML-based techniques. Once the ML model is trained, then the CPU's estimated value can be used by resource allocation schemes [10], which can dynamically assign the required resources to a particular VNF. To solve this resource requirement problem, we preset supervised ML models, Multiple linear regression (MLR), and Support Vector Machine (SVM) in comparison with Artificial Neural Network (ANN). All of the said models are implemented to solve an estimation problem, so they will be termed as multiple linear regression, support vector regression, and fitting neural network, respectively. Furthermore, they are evaluated against each other considering which model best fits the data represented by features.  
The main idea is to use a dataset for training our ML models, which can then estimate CPU value for each incoming traffic. For this purpose, first, the dataset must be constructed. Since our proposed models are based on supervised learning, therefore, the dataset must contain a pair of inputs (traffic vector) and outputs (real CPU utilization). The traffic vectors in dataset describes the incoming traffic to the VNFs, while the CPU utilization is the real-world CPU usage of these VNFs when the traffic is processed by them. The traffic is represented by 86 features (number of packets, IP address, total bytes, source and destination address, ports, etc., are a few mentions among others), which are extracted offline in 20 s batches (detail setup in Section 3). Once the dataset is constructed, it is further used to train our proposed ML models, which aim to learn the function that can accurately map the traffic features to CPU consumption. Figure 3 defines the considered problem in three layers.  
In this work, we modeled the behavior of real-world virtual network function that is shown on the bottom layer in Figure 3. It has been observed that the behavior of a particular VNF depends on different parameters, particularly its configuration, which, in our study, is considered to be constant. In system abstraction, which is the middle layer, the VNF is considered as block box, where the traffic represented by features enters the box and leaves the box, which in turn alters the VNF's resource requirements. The top layer represents machine learning architecture, where different regressors such as MLR, SVR, and FNN are used to estimate the CPU consumption of VNF. The considered ML models are only trained  
Figure 3: VNF Performance model.  
for CPU estimates of VNF because the exact configuration of VNF is hidden. The following section details the proposed models.
# Text Augmentation Using BERT for Image Captioning
## 3 Methodology
### Synonymous Augmentation  
Unlike image and speech processing augmentation by making random noises into the input signal (characters) is not suitable for text augmentation. The relative order of letters and their presence in a word can significantly affect the semantic meaning of the word itself. The best method of text augmentation is to rephrase sentence as a person would do. However, this approach is very complicated due to the size of a training datasets. One of the simpler, but still qualitative method, is synonymous replacing augmentation first introduced in [30].  
Let \(I\) be an image from a training set, \(C=\{c_{1},\dots,c_{k}\}\) be a set of captions corresponding to that image where each caption is a sequence of words \(c_{i}=(w_{i,1},w_{i,2},\dots,w_{i,l_{i}})\), \(l_{i}\) is a length of \(i\)-th caption. Additionally, fix some synonymous thesaurus \(T\) and let \(T(w_{i,j})=(s_{i,j,1},s_{i,j,2},\dots,s_{i,j,m_{i,j}})\) be a list of synonyms to a word \(w_{i,j}\) sorted in descending order of semantic closeness to the most frequently seen meaning of the word \(w_{i,j}\), \(m_{i,j}\)--the number of synonyms of the word \(w_{i,j}\) in the desaurus \(T\).  
The following operation is performed in order to generate a new caption \(c_{i}^{\prime}\) based on a \(c_{i}\). Fix some probability \(p\). For every word \(w_{i,j}\in c_{i}\), which has synonyms in a dictionary i.e., for which \(m_{i,j}>1\) replace it with synonym with a probability \(p\). To chose with which of the synonyms it will be replaced fix another probability \(q\). With a probability \(q\), replace a word with the most semantically close synonym \(s_{i,j,1}\). If the word was not replaced by the first synonym replace it with the second closest one \(s_{i,j,2}\) with a probability \(q\) and so on. Thus, the probability of replacing word with synonym \(s_{i,j,r}\) is equal to \(q^{r}\) and it exponentially decreases with the semantic similarity of the synonym to the original word. This operation of replacing a word with it synonym occurs independently for each word in a sentence.  
The above operation of obtaining a new caption based on an existing one is done \(d\) times, where \(d\) is called the augmentation coefficient. Accordingly, if the image has \(k\) captions after applying augmentation with a coefficient \(d\) it will have \(kd\) captions, i.e., training set is increased \(d\) times.
## Materials and Methods
### Two-Dimensional Fluid Model
#### 2.2.1 Model Equations
The rate of change of the electron energy density is described by [1]:  
\[\frac{\partial\omega_{e}}{\partial t}+\
abla\cdot\mathbf{\Gamma}_{\omega}+e \mathbf{\varepsilon}\cdot\mathbf{\Gamma}_{e}=R_{\omega} \tag{9}\]  
where \(\omega_{e}\) is the electron energy density, \(R_{\omega}\) is the energy loss or gain due to inelastic collisions, the term \(e\mathbf{\varepsilon}\cdot\mathbf{\Gamma}_{e}\) accounts for the ohmic or joule heating of the electrons in the electric field, and \(\mathbf{\Gamma}_{\omega}\) is the electron flux energy, that is described by:  
\[\mathbf{\Gamma}_{\omega}=\frac{5}{3}(-\mu_{e}\omega_{e}\mathbf{\varepsilon} -D_{\varepsilon}\
abla\omega_{e}) \tag{10}\]  
The electron energy loss or gain \(R_{\omega}\) is obtained by summing the collisional energy loss or gain over all reactions [14]:  
\[R_{\omega}=\sum_{j=1}^{P}x_{j}k_{j}N_{n}n_{e}\Delta\varepsilon_{j} \tag{11}\]  
where \(x_{j}\) is the mole fraction of the target species for reaction \(j\), \(k_{j}\) is the rate coefficient for reaction \(j\), \(N_{n}\) is the total neutral number density and \(\Delta\varepsilon_{j}\) is the energy loss from reaction \(j\). The electron energy density \(\omega_{e}\), the mean electron energy \(\varepsilon\), and the electron temperature \(T_{e}\) are correlated with each other through [30]:  
\[\omega_{e}=n_{e}\varepsilon=\frac{3}{2}k_{B}n_{e}T_{e} \tag{12}\]  
For non-electron species, the following equation was solved for the mass fraction of each species [30]:  
\[\rho\frac{\partial_{k_{B}}}{\partial t}+\rho(\boldsymbol{u}\cdot\
abla) \zeta_{k}=\
abla\cdot\boldsymbol{j}_{k}+R_{k} \tag{13}\]  
where \(\boldsymbol{j}_{k}\) is the diffusive flux vector, \(R_{k}\) is the rate expression for species \(k\), \(\boldsymbol{u}\) is the mass averaged fluid velocity vector, \(\rho\) denotes the density of the mixture and \(\zeta_{k}\) is the mass fraction of the \(k\)th species. The diffusive flux vector is defined as [30]:  
\[\boldsymbol{j}_{k}=\rho\zeta_{k}\boldsymbol{V}_{k} \tag{14}\]  
with \(\boldsymbol{V}_{k}\), being the multicomponent diffusion velocity for species \(k\). To initiate discharge in the reactor, electric potential should be applied between the electrodes, thus [PERSON]'s equation must also be considered in the model [14]:  
\[\varepsilon_{0}\
abla^{2}\varphi=-\sigma \tag{15}\]  
where \(\varphi\) is the electric potential, \(\varepsilon_{0}\) is the vacuum permittivity and \(\sigma\) is the charge density, that can be written in terms of density of the charged species, \(n_{k}\), and their charge, \(\epsilon\)Z\({}_{k}\)[31]:  
\[\sigma=e(\sum_{k=1}^{k}Z_{k}n_{k}-n_{e}) \tag{16}\]In this work, 16 different neutral and ionized species were considered in the model (Table 2). Thus, 16 continuity equations together with Poisson's equation were solved with the employment of a stabilized FEM.
# GNSS Timing Performance Assessment and Results Analysis
## 2 The Key Indicators and Assessment Method of GNSS Timing Performance
### Uncertainty Budget  
There are many kinds of errors during assessment of GNSS time transfer accuracy, including three parts: (1) the uncertainty from broadcast UTCO which is recorded as \(\sigma_{\text{UTCO}}\); (2) uncertainties related to satellite, signal propagation path, and receiver which are recorded as \(\sigma_{\text{GNSST}}\); (3) the uncertainty introduced by BIPM which is recorded as \(\sigma_{\text{BIPM}}\). These factors are not related; hence, the uncertainty components are independent. According to the principle of error analysis, the uncertainty of the results of GNSS time transfer accuracy assessment, recorded as \(\sigma\), is calculated as follows:  
\[\sigma=\sqrt{\sigma_{\text{GNSST}}^{2}+\sigma_{\text{UTCO}}^{2}+\sigma_{ \text{BIPM}}^{2}}. \tag{8}\]  
In Equation (7), \(\sigma_{\text{UTCO}}\) can be obtained from the provider's annual report, \(\sigma_{\text{BIPM}}\) can be obtained from BIPM, and \(\sigma_{\text{GNSST}}\) can be obtained by estimating the user equivalent range error (UERE) and time dilution of precision (TDOP), yielding the following formula [9; 10; 31].  
\[\sigma_{\text{GNSST}}=\text{UERE}\times\text{TODP}\div\text{C}, \tag{9}\]  
where C is the speed of light in m/s; TDOP can be replaced by the time transfer dilution of precision (TTDOP). TTDOP is related to the number of available satellites, about \(1/\sqrt{\text{N}}\).  
UERE is composed of two parts: the signal-in-space range error (SISE) which is associated with satellite clock error and orbit error, and the user equipment error (UEE) introduced by all residual error contributions in the range domain that are not under the direct control of the GNSS system. Considering the case that SISE and UEE are independent, UERE can be calculated by Equation (10), and the UERE evaluation is shown in Table 1[9; 10; 32; 33].  
\[\text{URERE}=\sqrt{\text{SISE}^{2}+\text{UEE}^{2}}. \tag{10}\]  
Because of the open environment in which the antenna of the monitoring receiver is mounted, the number of available satellites is 10 for BDS, eight for GPS, six for GLONASS, and six for Galileo. The uncertainty budget for the assessment results of GNSS time transfer accuracy is shown in Table 2. For GPS and GLONASS, UTCr-UTC(k) refers to UTCr-UTC(USNO) and UTCr-UTC(SU) respectively.  
It can be seen from Tables 1 and 2 that the uncertainty budget of GLONASS time transfer accuracy assessment is relatively large, which is due to the satellite clock error and orbit error. The \(\sigma_{\text{UTCO}}\) of BDS is large because the method of BDT traceability is indirect. First, the physical signal of BDT is controlled to the local physical reference. Then, BDT connects with UTC composed of two parts: one is the comparison link between the local physical reference and UTC(NTSC), while the other is the international network of time links between UTC(NTSC) and UTC. The uncertainty of BDT traceability is greater than that of other GNSSTs [28].
###### Abstract  
In this article, a novel statistical approach is proposed and applied to co-design signal and power integrity (SI/PI) in high-speed interconnects considering the non-linear power/ground noise generated by parallel buffers and bit-patterns. With increased data rates and decreased operating voltages, the allowed noise margin in high-speed interconnects is continuously reduced, and this trend requires SI/PI co-design. Specifically, non-linear power/ground noise associated with simultaneous switching circuits sharing a power delivery network (PDN) and bit-patterns must be carefully considered during the interconnects' design and analysis phase. In many cases, conventional electromagnetic (EM) and transient circuit simulators require heavy computational resources or even fail to deliver an accurate result. The proposed statistical method estimates the statistical eye-diagram in the high-speed interconnect considering power/ground noise and bit-patterns such as data bus inversion (DBI) coding. The accuracy and computational efficiency of the proposed method are validated by comparing the result with HSPICE transient simulation result. The proposed method is also compared with conventional statistical methods, such as peak distortion analysis (PDA) and statistical channel simulation in the transient simulator. Lastly, the proposed method is applied to the SI/PI co-design and co-analysis in the high bandwidth memory (HBM) interposer channel. Impacts of decoupling capacitors on hierarchical PDN impedance, statistical eye-diagram of the HBM channel, and bit error rate (BER) Bathtub curves are summarized. Finally, the BER eye-diagram is derived from the estimated statistical eye-diagram for timing and voltage analysis. The impacts of hierarchical PDN design and bit-patterns on SI/PI are discussed.  
statistical approach; signal/power integrity; high-speed interconnects
# Correction of Directional Effects in VEGETATION NDVI Time-Series
## 3 Methodology
### BRDF Correction
#### 3.2.3 Weights  
By assigning weights, we favor certain observations through the inversion process. The weights assigned to the ensemble of observations have both temporal and angular components. The analytical expression is given below:  
\[\sigma_{\text{TOC}}=0.5(c_{1}+c_{2}\,\rho_{\text{TOC}})\bigg{(}\frac{1}{\cos(1.058\,\theta_{\text{s}})}+\frac{1}{\cos(1.058\,\theta_{\text{p}})}\bigg{)} \tag{3}\]It was proposed by [31] to estimate uncertainties associated with TOC reflectances as a function of solar zenith (\(\theta_{s}\)) and viewing zenith (\(\theta_{v}\)) angles. The coefficients \(c_{1}\) and \(c_{2}\) in Equation (3) are wavelength-dependent and are taken from Table 1 in [31]. Thus, the angular weights are the inverse of the TOC uncertainty assessment with Equation (3). These weights allow setting less importance to observations that were acquired at large viewing and solar zenith angles.  
An effective temporal weight function is introduced by multiplying the prior covariance matrix (see Equation (10)) with a factor which is a function of the elapsed time since the last successful inversion. This temporal weight function has a decreasing exponential shape. It gives an enhanced weight to new observations. The timescale of such an exponential function is arbitrary set to 10 days, which means that a temporal weight obtained ten days before the last day of the accumulation period is 0.5. In contrast, the most recent observation has an effective temporal weight of 1 (maximum). This approach is consistent with the BRDF inversion methodology presented in [31; 32] and the Copernicus Climate Change Service's albedo product derived from VGT1 and VGT2 [[https://c3s.vgt.vito.be/documents/C3S_ALBEDO/D13.2-v1.0_ATBD_CDR_SA_SPOT-VGT_v1.0_PRODUCTS_v1.0.pdf](https://c3s.vgt.vito.be/documents/C3S_ALBEDO/D13.2-v1.0_ATBD_CDR_SA_SPOT-VGT_v1.0_PRODUCTS_v1.0.pdf)]([https://c3s.vgt.vito.be/documents/C3S_ALBEDO/D13.2-v1.0_ATBD_CDR_SA_SPOT-VGT_v1.0_PRODUCTS_v1.0.pdf](https://c3s.vgt.vito.be/documents/C3S_ALBEDO/D13.2-v1.0_ATBD_CDR_SA_SPOT-VGT_v1.0_PRODUCTS_v1.0.pdf)), accessed on 9 March 2021.
## 3 The OPTCON Algorithm
### The OPTCON3 Algorithm
 \tag{11}\]  
**Step I-2D:** Repeat steps [a]-[c] for each \(j=S+1,...,T\):[a): Calculate the components \(J_{D,j}\) and \(J_{C,j}\):  
\[J_{D,j} =\tfrac{1}{2}[x_{j}^{*\pi}-\tilde{x}_{j}]^{\prime}W^{xx}[x_{j}^{*\pi }-\tilde{x}_{j}]+[x_{j}^{*\pi}-\tilde{x}_{j}]^{\prime}W^{xx}[u_{j}^{*\pi}- \tilde{u}_{j}]\] \[\quad+\tfrac{1}{2}[u_{j}^{*\pi}-\tilde{u}_{j}]^{\prime}W^{uu}[u_{j}^{*\pi} -\tilde{u}_{j}]),\] (12) \[J_{C,j} =\tfrac{1}{2}tr(K_{j}^{xx}\Sigma^{\xi\xi}).\] [b): Calculate the matrix \(\Sigma^{\theta\theta}_{j/j}\):  
\[\Sigma^{\theta\theta}_{j/j}=\Sigma^{\theta\theta}_{j/j-1}-\Sigma^{\theta x}_{ j/j-1}(\Sigma^{xx}_{j/j-1})^{-1}\Sigma^{x\theta}_{j/j-1}\] (13) with \[\Sigma^{xx}_{j/j-1} =F^{x}_{\theta}\Sigma^{\theta\theta}_{j-1/j-1}(F^{x}_{\theta})^{ \prime}+\Sigma^{\varepsilon x}_{j}\] and \[\Sigma^{x\theta}_{j/j-1} =(\Sigma^{\theta x}_{j/j-1})^{\prime}=F^{x}_{\theta}\Sigma^{ \theta\theta}_{j-1/j-1},\quad\Sigma^{\theta\theta}_{j/j-1}=\Sigma^{\theta \theta}_{j-1/j-1}.\] [c): Compute the probing component \(J_{p}\):  
\[J_{P_{j}} =\tfrac{1}{2}tr[\Lambda^{xx}(\Lambda^{uu})^{-1}\Lambda^{ux} \Sigma^{xx}_{j/j}] \tag{15}\] \[\quad+tr[\Lambda^{xu}(\Lambda^{uu})^{-1}(B^{\prime}(K^{xx}D+K^{ \alpha\theta})+W^{ux}D)\Sigma^{\theta x}_{j/j}]\] \[\quad+\tfrac{1}{2}tr[((K^{\alpha\theta}+D^{\prime}K^{xx})B+D^{ \prime}W^{xx})(\Lambda^{uu})^{-1}\] \[\quad\times(B^{\prime}(K^{xx}D+K^{\alpha\theta})+W^{ux}D)\Sigma^{ \theta\theta}_{j/j}].\]  
**Step I-2E:** Calculate the sum of the deterministic, cautionary, and probing terms over the periods \(S,...,T\):  
\[\quad+\tfrac{1}{2}[u_{j}^{*\pi}-\tilde{u}_{j}]^{\prime}W^{uu}[u_{j}^{*\pi} -\tilde{u}_{j}]),\]0  
**Step I-2F:** Take a new control \(\overset{\circ}{u}_{S}=\overset{\circ\pi+1}{u}_{S}\) (a new point in the grid search) and go to step I-2A. **Step I-[3]:** Choose an optimal \(u^{*}_{S}\) with \(minJ=J^{*}_{d}(u^{*}_{S})\). End of grid search. Fix the corresponding \(x^{*}_{S}\). **Step II:** Calculate the following (a) and (b) for only one time period \(S\):  
(a) \[\quad+\tfrac{1}{2}[u_{j}^{*\pi}-\tilde{u}_{j}]^{\prime}W^{uu}[u_{j}^{*\pi} -\tilde{u}_{j}]),\]1 (16) and \[\quad+\tfrac{1}{2}[u_{j}^{*\pi}-\tilde{u}_{j}]^{\prime}W^{uu}[u_{j}^{*\pi} -\tilde{u}_{j}]),\]2 (b) \[\quad+\tfrac{1}{2}[u_{j}^{*\pi}-\tilde{u}_{j}]^{\prime}W^{uu}[u_{j}^{*\pi} -\tilde{u}_{j}]),\]3  
**Step III:** Update the parameter estimates \(\hat{\theta}\) and \(\Sigma^{\theta\theta}_{S/S}\):  
\[\quad+\tfrac{1}{2}[u_{j}^{*\pi}-\tilde{u}_{j}]^{\prime}W^{uu}[u_{j}^{*\pi} -\tilde{u}_{j}]),\]4  
\[\quad+\tfrac{1}{2}[u_{j}^{*\pi}-\tilde{u}_{j}]^{\prime}W^{uu}[u_{j}^{*\pi} -\tilde{u}_{j}]),\]5  
**Step IV:** Set \(\hat{\theta}=\hat{\theta}_{S/S}\) and \(\Sigma^{\theta\theta}=\Sigma^{\theta\theta}_{S/S}\), go to Step I and run the procedure for the time period \(S+1\). The OPTCON3 algorithm is finished when \(S=T\) and the approximately optimal dual control and state variables have been found for all periods.