Scientific Research

An Academic Publisher

**Non-Linear Tank Level Control for Industrial Applications** ()

Share and Cite:

*Applied Mathematics*,

**11**, 876-889. doi: 10.4236/am.2020.119057.

1. Introduction

1.1. Motivation

Tank level control is one of the most commonly encountered applications in industry [1] [2] [3]. A range of level control systems and methods are commonly used in industry. The control of liquid levels, for example in process tanks, is an important function. Examples include hot water storage tanks where water is periodically removed, and the level needs to be restored ready for the next wash cycle, effluent treatment plants, reaction vessels and steam boilers. The focus of this paper is on accurate liquid level control in single tank systems, where actuation (typically through a pump or proportional valve) can be modulated continuously (c.f. on/off type control actions). Specifically, accurate level control of plant where modulation of the level setpoint is also required is considered; such situations commonly arise in cascade control loops [4], for example, or situations in which supervisory Model Predictive Control (MPC) is implemented [5]. Such schemes can often occur in practice when a level control system is integrated into a wider plant Distributed Control System (DCS) employing supervisory optimization schemes above the plant control interface levels [1] [3] [5].

A complication of level control in such situation arises due to the non-linear equations of the level dynamics themselves. The typical industrial solution for regulatory control is to use a Proportional-Integral (PI) or Proportional Integral-Derivative controller which is tuned around a particular (fixed) operating point (tank level), using a linear approximation of the dynamics for tuning purposes. Although the linear approximation is accurate around the chosen operating point, if the system is perturbed away from this operating point (for example due to a disturbance or setpoint change), performance deteriorates rapidly. This arises since the linear approximation employed to tune the PI controller becomes less accurate due to non-linearity, and a fixed gain controller is not able to adapt or re-schedule its gains to match the new process conditions. The focus of this research is to develop a simple but effective solution to this issue.

1.2. Previous Work

A variety of advanced control schemes have been proposed to improve performance in such situations (e.g. non-linear MPC, adaptive/predictive control, gain scheduling, fuzzy logic control [5] - [10]). Although adaptive and adaptive/predictive controls such as [6] and [7] can be effective for tracking and adapting to slowly varying process parameters, a drawback lies in the time taken to adjust to quickly varying parameters; as this would be the case for level control after a setpoint change or disturbance, performance would deteriorate. An additional drawback of adaptive schemes is the added implementation complexity of parameter estimation and adaptation schemes, which may render industrial implementation (e.g. on a Programmable Logic Controller (PLC)) difficult (or impossible) in some cases (although recent work has proposed efficient schemes to overcome this, see e.g. [8] [9]). The computational overheads of non-linear MPC are also prohibitive for PLC implementation [5]. Gain scheduling schemes based around PI/PID controls have been found to perform well and have a much simplified implementation over adaptive schemes, but require a significant and time-consuming effort at the design stage to correctly tune. This occurs due to the need to linearize the non-linear tank dynamics at numerous operating points over the working range of the controller, with controller tuning to meet a specification at each point then required [10]. Fuzzy logic design often aims to mimic the operation of gain-scheduled approaches, and hence require similar effort at design time, but have heavier implementation overheads than gain schedulers [8] [9].

1.3. Contributions

In this paper, a very accurate technique based around feedback linearization and PI control is employed, in order to create a simple controller which can maintain linear performance over the full operating range of a uniform tank. Unlike previous gain-scheduled approaches, there is no significant effort required at design time, and unlike adaptive or fuzzy schemes, the implementation overhead is negligible. Taken together, it is argued that the approach is ideally suited for industrial implementation, for example on a PLC which may be operating within the framework of a DCS. The paper includes a discussion on the use of (linear) parameter estimation techniques within an adaptive control framework, although this clearly increases the implementation overheads. A simulation study coupled with experimental tests on a large-scale laboratory level control system using industrial control equipment illustrates the effectiveness of the proposed approach for both tracking and disturbance rejection.

1.4. Structure

The remainder of the paper is organized as follows. Section II introduces the system model, and provides an insight into the non-linearity of the process under consideration. Section III develops the non-linear control strategy, presents tuning formulae, discusses how linear parameter estimation may be employed within an adaptive framework and also describes implementation aspects. Section IV provides simulation-based and experimental validation of the approach, while Section V concludes the paper.

2. One Tank System Model

A schematic of the main components of a single uniform cross sectional tank is as depicted in Figure 1. For this tank, the following differential equation is easily obtained:

${A}_{T}\stackrel{\dot{}}{h}\left(t\right)={Q}_{in}\left(t\right)-{Q}_{out}\left(t\right)$ (1)

where A_{T} represents the cross-sectional area of the tank and
${Q}_{in}/{Q}_{out}$ represents the input and output flows respectively. Assuming that the input control signal is represented by
$u\left(t\right)$, the pump/valve has an effective flow gain
${K}_{v}$, and the output flow is through an aperture of effective area
${A}_{O}$, Equation (1) may be rewritten using Bernoulli’s equation as follows:

${A}_{T}\stackrel{\dot{}}{h}\left(t\right)={K}_{v}u\left(t\right)-{A}_{O}\sqrt{2gh\left(t\right)}$ (2)

Figure 1. Schematic of the single tank system.

Which can be re-written:

$\stackrel{\dot{}}{h}\left(t\right)=bu\left(t\right)-{a}^{\prime}\sqrt{h\left(t\right)}$ (3a)

$\text{With}:b=\frac{{K}_{v}}{{A}_{T}},\text{}{a}^{\prime}=\frac{{A}_{O}\sqrt{2g}}{{A}_{T}}$ (3b)

This is clearly a non-linear dynamic equation since the tank outflow depends upon the square root of the current tank level $h\left(t\right)$. However since $\sqrt{x}=x/\sqrt{x}$, a simple change of variables in (3) enables it to be put into the form of a linear equation, with a time-varying parameter $a\left(t\right)$ which is a function of the parameter ${a}^{\prime}$ and the current level $h\left(t\right)$ :

$\stackrel{\dot{}}{h}\left(t\right)=bu\left(t\right)-a\left(t\right)h\left(t\right)$ (4a)

$\text{With}:a\left(t\right)=\frac{{A}_{O}\sqrt{2g}}{{A}_{T}\sqrt{h\left(t\right)}}=\frac{{a}^{\prime}}{\sqrt{h\left(t\right)}}$ (4b)

Note that when
$h\left(t\right)=\text{1}$, the parameter
$a\left(t\right)$ is equivalent to
${a}^{\prime}$, and hence
${a}^{\prime}$ can be thought of as the nominal pole location for a unit level height. This pole location is effectively scaled by the reciprocal of the square root of the current level
$h\left(t\right)$. Denoting “s” as the usual time derivative operator, a 1^{st} order time-varying “transfer function” is easily created from (4):

$\frac{h\left(t\right)}{u\left(t\right)}=\frac{b}{s+a\left(t\right)}$ (5)

So it is seen that both the static gain and effective time constant of this 1^{st}-order process are functions of the instantaneous level
$h\left(t\right)$.

3. Non-Linear PI Control Strategy

3.1. 2-DOF PI Principle

Control of a first order system is easily achieved by the use of a PI controller. In order to maintain both tracking and regulation performance, to avoid the presence of an additional transmission zero in the forward path it is common to use a 2-DOF controller (see e.g. [11] for a general description and tuning rules for such controllers) such as that depicted in Figure 2.

This controller structure has the advantage of removing the controller zero from the closed loop dynamics if a regular error-activated PI controller was employed. It is advantageous in the current context as, with a simple modification, it produces linear closed-loop behavior for both tracking and regulation regardless of the current operating point $h\left(t\right)$. To see how this may be obtained, consider that for a first-order process, the closed loop dynamics when such a PI controller is employed is easily obtained through algebra and is given by:

$\frac{h\left(t\right)}{r\left(t\right)}=\frac{b{K}_{i}}{{s}^{2}+\left(a+b\right){K}_{p}s+b{K}_{i}}$ (6)

By setting $b{K}_{i}={\omega}_{n}^{2}$ and $\left(a+b\right){K}_{p}=2\zeta {\omega}_{n}$, one easily sees that (6) is simply the standard form of a unit-gain second-order transfer function. By comparing co-efficients, the controller is tuned for a desired damping ratio $\zeta $ and natural frequency ${\omega}_{n}$ as follows:

${K}_{p}=\frac{2\zeta {\omega}_{n}-a}{b}$ (7a)

${K}_{i}=\frac{{\omega}_{n}^{2}}{b}$ (7b)

3.2. Non-Linear Controller

Now, observe that only the model parameter b appears in (7b). Let us substitute for a our height-varying quantity $a\left(t\right)$ into (7a) and manipulate further:

${K}_{p}=\frac{2\zeta {\omega}_{n}-a\left(t\right)}{b}$ (8a)

$=\frac{2\zeta {\omega}_{n}}{b}-\frac{{a}^{\prime}}{b\sqrt{h\left(t\right)}}$ (8b)

Considering the required state feedback signal ${K}_{p}h\left(t\right)$ appearing at the summing junction to form $u\left(t\right)$ in Figure 2, and again using the identity $\sqrt{x}=x/\sqrt{x}$ this signal may be written as the combination of two sub-components:

Figure 2. 2-DOF Linear PI Controller.

${K}_{p}h\left(t\right)=h\left(t\right)\frac{2\zeta {\omega}_{n}}{b}-h\left(t\right)\frac{{a}^{\prime}}{b\sqrt{h\left(t\right)}}$ (9a)

$=h\left(t\right)\frac{2\zeta {\omega}_{n}}{b}-\sqrt{h\left(t\right)}\frac{{a}^{\prime}}{b}$ (9b)

Separation of these two components with opposing signs gives the simple non-linear PI controller structure as shown in Figure 3. The controller is effectively a regular PI controller-giving the control signal ${u}^{\prime}\left(t\right)$ as in Figure 2, with the addition of an extra linearizing control signal ${u}^{\u2033}\left(t\right)$ where the latter signal is formed as the feedback gain ${K}_{l}$ applied to the square root of the current level $h\left(t\right)$. This controller is tuned to give the desired damping ratio $\zeta $ and natural frequency ${\omega}_{n}$ for a (linear) second-order closed loop as follows:

${K}_{p}=\frac{2\zeta {\omega}_{n}}{b}$ (10)

${K}_{l}=\frac{{a}^{\prime}}{b}$ (11)

${K}_{i}=\frac{{\omega}_{n}^{2}}{b}$ (12)

In order to select appropriate values for natural frequency and damping ratio from natural design specifications, the following relations hold for determining these two parameters from settling time t_{s} (in seconds) and percent overshoot os (in %) for an underdamped second order system [1]:

$\zeta =\sqrt{\frac{\mathrm{ln}{\left(os/100\right)}^{2}}{{\pi}^{2}+\mathrm{ln}{\left(os/100\right)}^{2}}}$ (13)

${\omega}_{n}=\frac{4.6}{\zeta {\omega}_{n}}$ (14)

Closer observation of the control procedure reveals that it is a form of feedback linearization non-linear control [12], with a state-feedback controller with forward path integrator forming the control. Substituting the non-linear state feedback term ${u}^{\u2033}\left(t\right)={K}_{l}\sqrt{h\left(t\right)}$ into Equation (3) and using $u\left(t\right)={u}^{\prime}\left(t\right)\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{u}^{\u2033}\left(t\right)$ illustrates the method of loop linearization:

Figure 3. 2-DOF Non-linear PI Controller.

$\stackrel{\dot{}}{h}\left(t\right)=bu\left(t\right)-{a}^{\prime}\sqrt{h\left(t\right)}$ (15a)

$=b\left[{u}^{\prime}\left(t\right)+{u}^{\u2033}\left(t\right)\right]-{a}^{\prime}\sqrt{h\left(t\right)}$ (15b)

$=b\left[{u}^{\prime}\left(t\right)+{K}_{l}\sqrt{h\left(t\right)}\right]-{a}^{\prime}\sqrt{h\left(t\right)}$ (15c)

$=b{u}^{\prime}\left(t\right)+{a}^{\prime}\sqrt{h\left(t\right)}-{a}^{\prime}\sqrt{h\left(t\right)}$ (15d)

$=b{u}^{\prime}\left(t\right)$ (15e)

In other words, the non-linear state feedback cancels the effective pole and forces the process to behave as a weighted integrator of the control input term ${u}^{\prime}\left(t\right)$. Designing the linear state feedback term ${K}_{p}$ and integrator term ${K}_{i}$ is then trivial to achieve the required closed loop damping and natural frequency/settling time.

3.3. Self-Tuning/Adaptive Control Action

In some situations, it will be useful to implement a self-tuning and/or adaptive version of the non-linear controller. Adaptive controllers can be implemented using either direct or indirect approaches [13]; the latter is most commonly encountered. In an indirect approach, an estimation algorithm is first used to estimate the parameters of the controlled system on-line, and these parameters are used to calculate the controller gains. A block diagram of the typical layout of a self-tuning indirect adaptive controller is as shown in Figure 4. The adaptive controller in this context is formed by combining the on-line Exponentially Weighted Recursive Least Squares (EW-RLS) algorithm with the control law as described in the previous Section, using Equations (10), (11) and (12) to implement the tuning rules and map the parameter estimates to the controller gains, incorporating the design specification.

In the adaptive controller, the EW-RLS estimator is used to estimate the value of the process model parameters in real-time, and is appropriate for use in embedded adaptive control applications in which the model dimension is not excessive [13] [14]. In the general case, the EW-RLS algorithm is defined by the following expressions:

Figure 4. Indirect adaptive control scheme.

$\stackrel{^}{\beta}\left(t\right)=\stackrel{^}{\beta}\left(t-1\right)+K\left(t\right)e\left(t\right)$ (16)

$e\left(t\right)=y\left(t\right)-{x}^{\text{T}}\left(t\right)\stackrel{^}{\beta}\left(t-1\right)$ (17)

$K\left(t\right)=\frac{P\left(t-1\right)x\left(t\right)}{\lambda +{x}^{\text{T}}\left(t\right)P\left(t-1\right)x\left(t\right)}$ (18)

$P\left(t\right)=\frac{1}{\lambda}\left[P\left(t-1\right)-K\left(t\right){x}^{\text{T}}\left(t\right)P\left(t-1\right)\right]$ (19)

where $\beta \left(t\right)$ is the vector of estimated parameters of the system, $y\left(t\right)$ is the current measured output of the system under consideration, $x\left(t\right)$ is a vector of shifted previous input and output measurements of the system (regression variables), $K\left(t\right)$ is the estimator gain vector, $P\left(t\right)$ is the covariance matrix, $e\left(t\right)$ is the prior residual error and λ is the forgetting factor. In the context of this paper, only two parameters are required to be estimated, yielding a trivial implementation. From Equation (3), the two unknown parameters ${a}^{\prime}$ and b are linear in the measurements of $u\left(t\right)$ and $h\left(t\right)$, with the square of the latter being easy to calculate and its first derivative being obtainable numerically, e.g. using the trapezoidal rule [1]. Hence, although the dynamics are non-linear, the unknown parameters can be identified using EW-RLS directly.

3.4. Implementation

Equations (10), (11) and (12) above result in very simple design formulae to achieve the specified (linear) closed-loop behavior for the non-linear plant. In terms of a digital control implementation (e.g. on a PLC, microcontroller or within a DCS framework), several implementations are possible; the simplest is as follows. Applying Euler integration to the integral action of the controller, the following recursive equations result and can be employed for implementation:

$I\left(k\right)=I\left(k-1\right)+\frac{r\left(k\right)-h\left(k\right)}{{T}_{s}}$ (20)

$u\left(k\right)={k}_{i}I\left(k\right)-{k}_{p}h\left(k\right)+{k}_{l}\sqrt{h\left(k\right)}$ (21)

where integer $k$ represents discrete time index, $r\left(k\right)$ is the reference at step $k$, $h\left(k\right)$ is the measured height at step $k$, $I\left(k\right)$ is the integral of error at step $k$ and ${T}_{s}$ is the sample time (in seconds). The gains ${k}_{p}$, ${k}_{l}$ and ${k}_{i}$ are as given in Equations (10), (11) and (12). In terms of choice of sample time, given that the dynamics of most industrial level control applications are slow-moving, then multiples of one second will normally suffice, without causing undue computational burden. As with all physical control implementations, reset windup can be protected against by clamping the integrator output when the calculated control signal $u\left(k\right)$ reaches an upper or lower saturation limit. For details of the implementation of the EW-RLS algorithm, an efficient implementation for model structures such as Equation (3a) is given in reference [14].

4. Evaluation

In this section, a number of studies are described which were carried out to evaluate the effectiveness of the proposed non-linear controller both in simulation and on a large-scale experimental laboratory process rig.

4.1. Simulation Study

A simulation-based experiment was carried out on a small simulated tank within the Matlab®/Simulink® environment. The tank dynamics were assumed to be as follows:

$\stackrel{\dot{}}{h}\left(t\right)=2.0u\left(t\right)-0.5\sqrt{h\left(t\right)}$ (22)

The control design procedure as outlined in Section 3 was then applied to design a closed loop system with setting time of 13 seconds and 4.321% overshoot (natural frequency ${\omega}_{n}=0.\text{5}$ and optimal damping $\zeta =0.\text{7}0\text{71}$ ), to yield a closed-loop reference dynamics given by:

$\frac{0.25}{{s}^{2}+0.7071s+0.25}$ (23)

A Simulink model of the proposed controller and tank is as shown as Figure5. In the simulation, a series of three setpoints (50 mm, 100 mm and 150 mm) were successively applied at 0 seconds, 50 seconds and 100 seconds respectively. In addition, three small alternating disturbance flows were also introduced at 25 seconds, 75 seconds and 125 seconds. The total simulation time was 150 seconds. The results obtained are as shown in Figure6. The Figureshows the commanded setpoint level (black trace) versus time, and also the simulated actual level (red trace) also versus time. For this Figure, it can clearly be observed that linear behaviors are achieved across the range of commanded levels, with the required overshoot and setline time being achieved in each case. Disturbance rejection is quick, effective and uniform at each of the three operating levels.

Figure 5. Matlab®/Simulink® model of the controller.

Figure 6. Simulated tracking and disturbance rejection.

4.2. Experimental Study

In order to further evaluate the proposed control method described in the previous Section on real industrial hardware, a second experiment was carried out on an existing test facility at Teesside University, as shown in Figure 7. The configuration of the rig for the purposes of this paper was to control the level of water in the holding tank using an inverter-connected pump, with the drain outflow value to the sump set at 50% open. Level was instrumented using an accurate float sensor. Control signals to and from the rig were implemented by 4 - 20 mA current loops to a local Schneider® PLC, with TIA-232 link to a local PC executing Simulink© for data acquisition purposes. The Modbus protocol allows setpoint, control signal and level sensor signals to be read/written in real-time from a Simulink® model executing within the Real-Time Windows Target environment. For generation of control code for embedded targets directly from Simulink© models, the interested reader is referred to [15], for example.

The control cabinet and PLC are depicted in Figure 8, which also shows the inverter driving the pump motor, power supply/protection and analog interfaces. Previous experiments on System Identification have been carried out on this rig, and w.r.t. Equation (3), the following model parameters have been identified from a drain-and-fill test:

$\stackrel{\dot{}}{h}\left(t\right)=0.0122u\left(t\right)-0.00497\sqrt{h\left(t\right)}$ (24)

Comparison of the model tank level with the real tank level during the drain-and-fill test is as shown in Figure 9. The calibrated model fits the data with accuracy > 90%, with the main deviation occurring with the tank being close to empty, away from the main operating points of the test rig.

The control design procedure as outlined in Section 3 was then applied to design a closed loop system with setting time 100 seconds and 4.321% overshoot (natural frequency ${\omega}_{n}=0.0\text{65}$ and optimal damping $\zeta =0.\text{7}0\text{71}$ ), to yield a closed-loop reference dynamics given by:

Figure 7. Experimental test rig.

Figure 8. Experimental test rig control cabinet.

Figure 9. Model level comparison to real level.

$\frac{0.25}{{s}^{2}+0.7071s+0.25}$ (25)

The settling time was chosen to be 100 seconds to prevent problems due to saturation of the pump, as values significantly shorter than this lead to saturation and integral action clamping. The controller was then implemented within the PLC, using a sample time of one second. A function block implementation was created, implementing Equations (16) and (17). Only a very small amount of instructions (analog read/write operations, addition, multiplication, square root and memory read/write operations) were required. An experiment was carried out in which two-level setpoints (0.5 m and 1.0 m) were applied at 0 seconds and 250 seconds. The total simulation time was 500 seconds. The results obtained are as shown in Figure 10.

The Figureshows the commanded setpoint level (black trace) versus time, and also the simulated actual level (red trace) also versus time. For this Figure, it can clearly be observed that linear behavior was achieved across the two commanded levels, with the required overshoot and settling time being (approximately) achieved in each case, despite the presence of some noise affecting the level measurement.

5. Summary

The focus of this paper has been upon accurate liquid level control in single tank systems which can be actuated continuously, and modulation of the level setpoint is also required. A simple but accurate technique based around feedback linearization and Proportional Integral (PI) control has been introduced. It has been shown that this simple controller can maintain linear performance over the full operating range of a uniform tank, using only parameters estimated from a simple fill-and-drain experimental test. Implementation details have been discussed. Simulation results and experimental results on a large-scale laboratory level control system implemented using an industrial PLC illustrate the effectiveness of the proposed approach for both level tracking and disturbance rejection.

Figure 10. Experimental tracking of level setpoint.

It is concluded that the proposed scheme can potentially benefit industrial level control applications, both for stand-alone operations and those implemented within a wider supervisory control scheme within a DCS. Future work will focus on the application of feedback linearization to multiple chained tank applications.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

[1] | Dorf, R.C. and Bishop, R.H. (2004) Modern Control Systems. 10th Edition, Prentice-Hall, Englewood Cliffs, NJ. |

[2] |
Smith, R.S. and Doyle, J. (1988) The Two Tank Experiment: A Benchmark Control Problem. Proceedings of the American Control Conference, Atlanta, 15-17 June 1988, 403-415. https://doi.org/10.23919/ACC.1988.4790058 |

[3] |
Liptak, B.G. (2005) Instrument Engineers’ Handbook Volume Two: Process Control and Optimization. 4th Edition, CRC Press, Boca Raton, 2464 p. https://doi.org/10.1201/9781420064001 |

[4] |
Bolton. W. (2015) Instrumentation and Control Systems. 2nd Edition, Elsevier, Amsterdam. https://doi.org/10.1016/B978-0-08-100613-9.00004-3 |

[5] | Camacho, E.F. and Bordons, C. (2004) Model Predictive Control. 2nd Edition, Springer Verlag, Berlin. |

[6] |
Cartes, D. and Wu, L. (2005) Experimental Evaluation of Adaptive Three-Tank Level Control. ISA Transactions, 44, 283-293. https://doi.org/10.1016/S0019-0578(07)60181-5 |

[7] |
Short, M. and Abugchem, F. (2017) A Microcontroller-Based Adaptive Model Predictive Control Platform for Process Control Applications. Electronics, 6, 88. https://doi.org/10.3390/electronics6040088 |

[8] | Xiao, Q.H., Zou, D.Q. and Wei, P. (2010) Fuzzy Adaptive PID Control of Tank Level. 2010 International Conference on Multimedia Communications, Hong Kong, 7-8 August 2010, 149-152. |

[9] | Sislin, R., Da Silva, F.V., Gedrite, R., Jokinen, H. and Rajan, D.K. (2016) Mathematical Modeling and Development of a Low Cost Fuzzy Gain Schedule Neutralization Control System. Engineering Letters, 26, 353-357. |

[10] |
Wang, M. and Crusca, F. (2002) Design and Implementation of a Gain Scheduling Controllerfor s Water Level Control System. ISA Transactions, 41, 323-331. https://doi.org/10.1016/S0019-0578(07)60091-3 |

[11] | Vilanova, R., Alfaro, V.M. and Arrieta, O. (2011) Analytical Robust Tuning Approach for Two-Degree-of-Freedom PI/PID Controllers. Engineering Letters, 19, 204-214. |

[12] | Slotine, J.J. and Li, W. (1988) Applied Non-Linear Control. Pearson Education, London. |

[13] | Astrom, K.J. and Wittenmark, B. (1995) Adaptive Control. 2nd Edition. Addison Wesley, Boston. |

[14] |
Short, M. (2012) Fast Online Identification of Low-Order Time-Delayed Industrial Processes. Electronics Letters, 48, 152-153. https://doi.org/10.1049/el.2011.3400 |

[15] | Kamiyama, T., Tamura, M., Soeda, T., Yoo, M. and Yokoyama, T. (2012) An Embedded Control Software Development Environment with Simulink Models and UML Models. IAENG International Journal of Computer Science, 39, 261-268. |

Copyright © 2020 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.