\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[english]{babel}



\usepackage{float}






%%% Your header encountered errors when converting with LaTeXML
%%% Please uncomment it incrementally and check for problems
%%%   If you need assistance, contact us at help@authorea.com
%%% 
%
%\iflatexml
%\documentclass[11pt, letterpaper]{article}
%\fi
%\usepackage{graphicx}
%
%\usepackage{amsmath}
%
%\usepackage{amssymb}
%
%\usepackage{setspace}
%
%% \usepackage{epstopdf} % Not (yet) supported by LaTeXML.
%
%% \usepackage{empheq} % Not (yet) supported by LaTeXML.
%
%\usepackage[letterpaper,left=1in,right=1in,top=1in,bottom=1in]{geometry}
%
%% \usepackage[definethebibliography]{easybib} % Not (yet) supported by LaTeXML.
%
%\citestyle{C}
%\usepackage{float}
%
% \usepackage{color}
%\renewcommand{\refname}[][]{Literature Cited}
%\newtheorem{remark}{Remark}
%\newtheorem{lemma}{Lemma}
%\newtheorem{proposition}{Proposition}
%\newtheorem{corollary}{Corollary}
%\newtheorem{theorem}{Theorem}
%\newtheorem{alg}{Algorithm}
%\usepackage{url}
%\newcommand{\m}[][]{\mathbb}
%\allowdisplaybreaks
%\iflatexml
%\documentclass[11pt, letterpaper]{article}
%\fi
%\usepackage{graphicx}
%\usepackage{amsmath}
%\usepackage{amssymb}
%\usepackage{setspace}
%% \usepackage{epstopdf} % Not (yet) supported by LaTeXML.
%% \usepackage{empheq} % Not (yet) supported by LaTeXML.
%\usepackage[letterpaper,left=1in,right=1in,top=1in,bottom=1in]{geometry}
%% \usepackage[definethebibliography]{easybib} % Not (yet) supported by LaTeXML.
%\usepackage{float}
%\usepackage{color}
%\renewcommand{\refname}[][]{Literature Cited}
%\usepackage{url}
%\newcommand{\m}[][]{\mathbb}
%
%\setstretch{1.5}
%\newpage
%\section*{Introduction}
%The pulp and paper industry has a profound influence on the economy of the world, which produces pulp, paper, paperboard, and various cellulose-based products.
%Even though electronic media and paperless communication have been widely expanded, the global pulp and paper market is growing steadily at a rate of over 1\% per year (missing citation); (missing citation). For instance, the global consumption of paper and board amounted to an estimated 399 million metric tons in 2020. It is expected that demand will increase steadily over the next decade, reaching approximately 461 million metric tons in 2030 (missing citation).
%The vast majority of increasing demand mainly comes from the following categories of products (missing citation): personal hygiene paper products, food packaging products, corrugated packaging materials, and paper-based medicinal materials (shown in Figure \ref{paper}), which closely related to the growing e-commerce business and awareness of safety and hygiene, especially under the pandemic situation of COVID-19. 
%Advanced process control and state-of-the-art process optimization techniques would provide enormous opportunities for maximized efficiency and optimized energy usage to satisfy the steadily growing need.
%The pulp and paper mills aim to convert wood chips into pulp, the raw material for different types of products. In general, a pulping process can be classified into mechanical pulping, chemical pulping, and thermomechanical pulping according to the fundamental mechanism of separating wood fibers (missing citation). More than two-thirds of globally produced pulp comes from Kraft pulping which belongs to the chemical pulping process (missing citation). In the Kraft pulping process, the conversion of wood chips into pulp mainly takes place in a pressured vertical cylindrical reactor known as the pulp digester, which operates in a batch manner or as a continuous process. Due to lower space requirements and lower energy costs, continuous pulp digesters are predominantly used in industrial applications (missing citation). As illustrated in Figure \ref{digester}, the typical continuous digesters consist of three basic zones: an impregnation zone, a cook zone, and a wash zone. White liquor and pre-steamed wood chips are introduced at the top of the digester into the impregnation zone where the liquor penetrates the wet chips. After that, the two streams flow downwards into the cook zone where the most delignification reactions occur. Then, the spent liquor is withdrawn from the digester at extraction screens. At the same time, the cooked pulp moves downwards to the wash zone where the chips are washed by the counter-current flow of cold wash liquor. The cooked pulp is then removed from the bottom of the digester. In particular, the cooking degree is evaluated by the Kappa number, which denotes the residual lignin content of the pulp.
%Due to the complex nature of the delignification process, numerous mathematical models for the pulp digesters have been proposed over the last few decades. Three widely used dynamic models are known as the Purdue, Gustafson, and Andersson models, which have similar conceptual bases and assumptions (missing citation); (missing citation); (missing citation). 
%These kinetic models show the effect of a change of temperature on the reaction rate constants for different wood components according to the Arrhenius expression (missing citation). The main differences are the numbers of wood components, the assumption of consecutive or parallel reactions, and the assumption about how the delignification reaction takes place along the digester length. Among them, the most commonly used model is the extended Purdue model, which is followed and further developed by other researchers (missing citation); (missing citation); (missing citation).\\
%When it comes to the controller design for the pulp digesters, a wide variety of control methods, such as reduced dimension control, robust control, and generic algorithms have been investigated  (missing citation); (missing citation). Since one of the main objectives of the pulp digester is to produce pulp with specific properties with minimum chemical and energy inputs, the digester is usually operated in a constrained setting(missing citation). 
% In this case, model predictive control (MPC) as a widely deployed methodology in different types of  processes is capable of handling such requirements.  Michaelsen et al. (missing citation) developed a model predictive controller for a Kamyr digester using a real-time mechanistic model compensated by an optimal state estimator and the performance demonstrated to be superior when compared with proportional-integral (PI) control in offline simulations. 
%Lee and Datta  (missing citation) designed a model predictive control  system coupled with an extended Kalman filter for a batch-type pulp digester.
%Along the same line of work, Wisnewski and Doyle (missing citation) analyzed the performance of linear MPC and nonlinear MPC for set-point tracking and unmeasured disturbance rejection.
%Alexandridis and Sarimveis (missing citation) employed an adaptive MPC based on a radial basis function ANN model for Kappa number control of continuous pulp digesters.
%H.-K. Choi and J. S.-I. Kwon (missing citation); (missing citation); (missing citation); (missing citation) developed a class of MPC for continuous pulp digester and batch pulp digester based on the proposed multi-scale model.
%However, even though the aforementioned works have made a valuable contribution toward the modelling and controller design for the pulp digester, there still some aspects which did not receive much attention and/or accurate consideration. First of all, most of the works depend on the approximation of the PDE into a large-scale ordinary differential equation (ODE), which is generally prone to approximation error and is difficult to capture the spatial kinetics properties within the digester. In addition, these approaches need to perform the spatial discretization in an early lumping manner, which might induce numerical instability and/or alter the fundamental control theoretical properties, such as controllability, observability, and stabilizability (missing citation); (missing citation). Secondly, the target of the kappa number and cooking temperature setting are considerably different between softwood and hardwood pulp, the shift operations, therefore, need to be taken into account when the controller is developed. Thirdly, the co-current flow and counter-current flow in the digester need to be considered in accurate manner, mainly due to the fact that the main control manipulation is performed in the top part of the digester, while the controlled variable is physically measured at the outlet of the digester (missing citation).
%Motivated by the aforementioned issues, this work considers a bilinear transformation of a continuous infinite-dimensional system to a discrete one by the application of the Cayley-Tustin time discretization approach, in which the physical characteristics (energy) and theoretical properties of considered systems are preserved (missing citation). In addition, the finite-dimensional MPC setting is extended to the distributed parameter pulp digester system to realize target tracking and achieve optimal grade transitions. The cascade zones (cook zone and wash zone) of the digester are considered as an extended distributed parameter system, which is described by 10 hyperbolic PDEs.
%In this manuscript, the claimed novelty is the extension of linear MPC designs for the finite-dimensional system to the infinite-dimensional one to realize target tracking. Particularly, the pulp digester system described by a set of coupled nonlinear hyperbolic PDEs is considered, which can capture the spatiotemporal evolution of wood chips and wash liquor. In addition, the linearized continuous-time model is transformed to a linear discrete-time infinite-dimensional model, and the fundamental continuous-time properties (including stability, controllability and observability) are preserved under the Cayley–Tustin transformation. Moreover, a discrete-time observer for the system of hyperbolic PDEs is proposed, accounting for the available output measurement and the system states reconstruction. Finally, the controller design provides optimal asymptotic stabilization and target tracking realization of the system with the inclusion of output and input constraints.
%The article is organized as follows: we initially present a dynamic model to describe the delignification process in the continuous pulp digester. Then, the developed model is linearized and discretized in time by utilizing the Cayley-Tustin approach. Lastly, we propose a Luenberger observer-based model predictive controller for the discrete model to realize the target tracking. The simulation results demonstrate the performance of the proposed controller design.
%\section*{Problem Formulation}
%In this section, the mathematical model of the pulp digester process is introduced. In particular, the equilibrium profiles are calculated to proceed with the linearization of the original nonlinear system. For the sake of simplicity, an infinite-dimensional representation is given to illustrate the linearized system. Finally, a discrete-time model is obtained utilizing the Cayley-Tustin discretization framework.
%\subsection*{Model description}
%The main section of a continuous pulp digester can be divided into the co-current zone and counter-current zone, which are also referred to as cook zone and wash zone respectively, as shown in Figure \ref{scheme}. Each zone can be modelled by a set of nonlinear coupled PDEs, which are derived from the first principle. The models of cook zone and wash zone are similar, but different flow direction of the liquor. Generally, each volume in the digester contains three phases, solid phase, entrapped liquor phase, and free liquor phase. In particular, the combination of solid phase and entrapped liquor phase is further referred to as the wood chip phase, and the details of the conceptual model of the mass in the digester are shown in Figure \ref{phase}. The wood chips are often assumed to consist of five components, and the entrapped liquor and free liquor contain four same components. For more details of the model, one can see the references(missing citation); (missing citation).
%Based on the model proposed by Michelsen, the following assumptions are considered to obtain a simple model which describes some dominant dynamics of the process. Wood chips are assumed to be of equal size and to have a constant volume during cooking. The volume flows of wood chips and free liquor are assumed to be equal at all space positions, and the chip compaction profile is assumed to be stationary, expressed as a piecewise linear function of position. The cross-sectional area is assumed to be constant in the digester. For simplicity, the degradation of carbohydrates is assumed to follow a lignin reduction linearly and the concentrations of dissolved solids in the entrapped liquor and free liquor can be omitted  (missing citation); (missing citation); (missing citation). According to Gustafson et al. (missing citation), the carbohydrate reaction rate in the bulk delignification period is proportional to the lignin reaction rate.
%Hence, from control and monitoring point of view, the following variables are important, which are considered as the states of the model:
%$\rho_{s} (z, t)$: concentration of lignin of solid phase (\%)
%$\rho_{e} (z, t)$: concentration of alkali of entrapped liquor phase ($g/l$)
%$\rho_{f} (z, t)$: concentration of alkali of free liquor phase ($g/l$)
%$T_c (z, t)$: temperature of wood chip phase (K)
%$T_f (z, t)$: temperature of free liquor phase (K)
%Based on the abovementioned assumption, the mathematical models of cook zone and wash zone can be simplified to the following equations:
%\begin{equation}
%\begin{aligned}
%&0\leq z\leq L_{co}\left\{
%\begin{array}
%{l}
%\frac{\partial \rho_{s}^{\mathcal{C}} (z, t)}{\partial t}=-v_c\frac{\partial \rho_{s}^{\mathcal{C}}  (z, t)}{\partial z}-R_{s}^{\mathcal{C}}\\ 
%\frac{\partial \rho_{e}^{\mathcal{C}} (z, t)}{\partial t}=-v_c\frac{\partial  \rho_{e} ^{\mathcal{C}} (z, t)}{\partial z}+\frac{D_{EA}^{\mathcal{C}}}{\varepsilon _c}(\rho_{f}^{\mathcal{C}} (z, t)-\rho_{e}^{\mathcal{C}} (z, t))-\rho_{ODW}(b_1+0.47b_2)R_{s}^{\mathcal{C}}\\ 
%\frac{\partial \rho_{f}^{\mathcal{C}}  (z, t)}{\partial t}=-v_f\frac{\partial \rho_{f}^{\mathcal{C}}  (z, t)}{\partial z}-\frac{D_{EA}^{\mathcal{C}}}{\varepsilon _f}(\rho_{f}^{\mathcal{C}} (z, t)-\rho_{e}^{\mathcal{C}} (z, t))\\ 
%\frac{\partial T_c^{\mathcal{C}} (z, t)}{\partial t}=-v_c\frac{\partial T_c^{\mathcal{C}} (z, t)}{\partial z}+\left [ 1.47\rho_c \Delta H_R R_{s}^{\mathcal{C}}+U(T_f^{\mathcal{C}} (z, t)-T_c^{\mathcal{C}} (z, t))\right ] / C_{pe} \\
%\frac{\partial T_f^{\mathcal{C}} (z, t)}{\partial t}=-v_f\frac{\partial T_f^{\mathcal{C}} (z, t)}{\partial z}-U(T_f^{\mathcal{C}} (z, t)-T_c^{\mathcal{C}} (z, t)) /C_{pf}+b(z)u(t)
%\end{array}
%\right.\\
%&0\leq z \leq L_{ws}\left\{
%\begin{array}
%{l}
%\frac{\partial \rho_{s} ^{\mathcal{W}}(z, t)}{\partial t}=-v_c\frac{\partial \rho_{s} ^{\mathcal{W}}(z, t)}{\partial z}-R_{s}^{\mathcal{W}}\\ 
%\frac{\partial \rho_{e}^{\mathcal{W}}(z, t)}{\partial t}=-v_c\frac{\partial  \rho_{e}^{\mathcal{W}} (z, t)}{\partial z}+\frac{D_{EA}^{\mathcal{W}}}{\varepsilon _c}(\rho_{f}^{\mathcal{W}}(z, t)-\rho_{e}^{\mathcal{W}}(z, t))-\rho_{ODW}(b_1+0.47b_2)R_{s}^{\mathcal{W}}\\ 
%\frac{\partial \rho_{f} ^{\mathcal{W}}(z, t)}{\partial t}=v_f\frac{\partial \rho_{f}^{\mathcal{W}}(z, t)}{\partial z}-\frac{D_{EA}^{\mathcal{W}}}{\varepsilon _f}(\rho_{f}^{\mathcal{W}}(z, t)-\rho_{e}^{\mathcal{W}}(z, t))\\ 
%\frac{\partial T_c^{\mathcal{W}}(z, t)}{\partial t}=-v_c\frac{\partial T_c^{\mathcal{W}}(z, t)}{\partial z}+\left [  1.47\rho_c \Delta H_R R_{s}^{\mathcal{W}}+U(T_f^{\mathcal{W}}(z, t)-T_c^{\mathcal{W}}(z, t))\right ] / C_{pe} \\
%\frac{\partial T_f^{\mathcal{W}}(z, t)}{\partial t}=v_f\frac{\partial T_f^{\mathcal{W}}(z, t)}{\partial z}-U(T_f^{\mathcal{W}}(z, t)-T_c^{\mathcal{W}}(z, t)) /C_{pf}
%\end{array}
%\right.
%\label{nonpde} 
%
%\end{aligned}
%
%\end{equation}
%where the superscript $\mathcal{C}$ and $\mathcal{W}$ indicate the cook zone and wash zone, respectively. $R_{s}^{\mathcal{C}}$ and $R_{s}^{\mathcal{W}}$ are the rate of consumption of mass of solid per chip volume in cook zone and wash zone, describing by the Arrhenius coefficients, $R_{s}^{\mathcal{C}}=A_{1}e^{-\frac{E_{1}}{T_c^{\mathcal{C}}}}\rho_{e}^{\mathcal{C}}(\rho_{s}^{\mathcal{C}}-\rho_{s}^{0 })$, $R_{s}^{\mathcal{W}}=A_{1}e^{-\frac{E_{1}}{T_c^{\mathcal{W}}}}\rho_{e}^{\mathcal{W}}(\rho_{s}^{\mathcal{W}}-\rho_{s}^{0 })$. $\rho_{s}^{0 }$ denotes the non-reactive lignin in wood. The mass diffusion rate $D_{EA}$ usually take the form, $D_{EA}^{\mathcal{C}}=0.16\sqrt{T_c^{\mathcal{C}}}e^{\frac{-2452.4}{T_c^{\mathcal{C}}}}(-2.0\rho_{s}^{\mathcal{C}}+0.13(\frac{\rho_{e}^{\mathcal{C}}}{32.0})^{0.55}+0.58)$ and $D_{EA}^{\mathcal{W}}=0.16\sqrt{T_c^{\mathcal{W}}}e^{\frac{-2452.4}{T_c^{\mathcal{W}}}}(-2.0\rho_{s}^{\mathcal{W}}+0.13(\frac{\rho_{e}^{\mathcal{W}}}{32.0})^{0.55}+0.58)$. $\varepsilon_c$ is the chip compaction, which increases from the entry through the cook zone, reaching a maximum at the main extraction, $\varepsilon_c(z)=\varepsilon _{10}+\varepsilon _{11}z$, and $\varepsilon_f(z)=1-\varepsilon_c(z)$. $U$ denotes the heat-transfer coefficient, and the heat capacities of the entrapped and free liquor phases ($C_{pe}$ and $C_{pf}$) are determined by the mixing rules based on weighted averages (missing citation). 
%In this case, the input is the temperature of free liquor and the distribution is described by the actuation distribution function $b(z)=1_{[z_a-\varepsilon, z_a+\varepsilon]}$. 
%The controlled output is defined as the concentration of lignin of solid phase at the outlet in wash zone, and the measurements of each state at the outlet are assumed available which are taken as the measured outputs.
%The boundary conditions for the cook zone are given at $z=0$ (see Figure \ref{scheme}):
%\begin{align}
%\rho_{s}^{\mathcal{C}} (0,t)&=\rho_{s_0}, \rho_{e}^{\mathcal{C}} (0,t)=\rho_{e_0}, \rho_{f}^{\mathcal{C}} (0,t)=\rho_{f_0}, 
%T_c^{\mathcal{C}} (0, t)=T_{c0}, T_f^{\mathcal{C}} (0, t)=T_{f0}
%
%\end{align}
%and for the wash zone at  $z=0$ and $z=L_{ws}$:
%\begin{small}
%\begin{align}
%\rho_{s}^{\mathcal{W}}(0,t)&=\rho_{s}^{\mathcal{C}}(L_{co},t),  \rho_{e}^{\mathcal{W}}(0,t)=\rho_{e}^{\mathcal{C}}(L_{co},t),
%\rho_{f}^{\mathcal{W}}(L_{ws}, t)=\rho_{f_0}, 
%T_c^{\mathcal{W}}(0, t)&=T_c^{\mathcal{C}}(L_{co}, t),T_f^{wash}(L_{ws}, t)=T_{f0}
%\label{bcwash} 
%
%\end{align}
%
%\end{small}
%\subsection*{Model Linearization}
%In order to reduce the mathematical complexity of the nonlinear process model, we make the linearization around the system equilibrium point or spatially nonuniform steady-state by setting the time derivative terms to be zero. The steady-states are solved numerically using the finite difference method. As a result, the corresponding steady-states profiles are illustrated in Figure \ref{ss}.
% Additionally, we assume that $L_{co}$ and $L_{ws}$ are the same, which denote as $L$ for notation simplicity, and introduce the following notations to define the local dynamics from the states of the cook zone to the states of the wash zone:
%\begin{equation}
%\begin{aligned}
%x^{(\mathcal{C})}(z,t)&=[\rho_{s}^{\mathcal{C}}(z,t); \rho_{e}^{\mathcal{C}}(z,t); \rho_{f}^{\mathcal{C}}(z,t); T_c^{\mathcal{C}}(z,t); T_f^{\mathcal{C}}(z,t)]\\
%x^{(\mathcal{W})}(z,t)&=[\rho_{s}^{\mathcal{W}}(z,t); \rho_{e}^{\mathcal{W}}(z,t); \rho_{f}^{\mathcal{W}}(z,t); T_c^{\mathcal{W}}(z,t); T_f^{\mathcal{W}}(z,t)]
% 
%\end{aligned}
%
%\end{equation}
%The above nonlinear model Equation (\ref{nonpde}) can be linearized applying the Taylor series expansion. The linearized model is then obtained as:
%\begin{equation}
%\begin{aligned}
%\dot{x}^{(\mathcal{C})}(z,t)&=\mathcal{A}_1x^{(\mathcal{C})}(z,t)+\mathcal{B}_1u(t)\\
%\dot{x}^{(\mathcal{W})}(z,t)&=\mathcal{A}_2x^{(\mathcal{W})}(z,t)\\
%y_c(t)&=\mathcal{C}_2x^{(\mathcal{W})}(z,t)\\
%y_{m}(t)&=
%\begin{bmatrix}
%\mathcal{C}_{m1}x^{(\mathcal{C})}(z,t)\\
%\mathcal{C}_{m2}x^{(\mathcal{W})}(z,t)
%
%\end{bmatrix}
%\label{sspace} 
%
%\end{aligned}
%
%\end{equation}
%where $z_1\in [0,L]$ and $z_2\in [0,L]$. The operator $\mathcal{A}_1(\cdot)=(V_1\frac{\partial }{\partial z}+\psi_1(z))(\cdot)$ and $\mathcal{A}_2(\cdot)=(V_2\frac{\partial }{\partial z}+\psi_2(z))(\cdot)$, where $V_1=\text{diag}\{-v_c, -v_c, -v_f, -v_c, -v_f\}$ and $V_2=\text{diag}\{-v_c, -v_c, v_f, -v_c, v_f\}$ and
%\begin{small}
%\begin{align}
%\psi_1(z)&=
%\begin{bmatrix}
%J_{11}(z) &  J_{12}(z) & 0 & J_{14}(z) &0 \\ 
%J_{21}(z) & J_{22}(z) & J_{23}(z) & J_{24}(z) & 0  \\ 
% J_{31}(z)&  J_{32}(z) & J_{33}(z) & J_{34}(z) &  0 \\ 
%  J_{41}(z)&  J_{42}(z) & 0 & J_{44}(z) & J_{45}(z) \\ 
%0& 0 & J_{53}(z) & J_{54}(z) & J_{55}(z) 
%\end{bmatrix}
%, \psi_2(z)=
%\begin{bmatrix}
%\bar{J}_{11}(z) &  \bar{J}_{12}(z) & 0 & \bar{J}_{14}(z) &0 \\ 
%\bar{J}_{21}(z) & \bar{J}_{22}(z) & \bar{J}_{23}(z) & \bar{J}_{24}(z) & 0  \\ 
% \bar{J}_{31}(z)&  \bar{J}_{32}(z) & \bar{J}_{33}(z) & \bar{J}_{34}(z) &  0 \\ 
%  \bar{J}_{41}(z)&  \bar{J}_{42}(z) & 0 & \bar{J}_{44}(z) & \bar{J}_{45}(z) \\ 
%0& 0 & \bar{J}_{53}(z) & \bar{J}_{54}(z) & \bar{J}_{55}(z) 
%\end{bmatrix}
%\nonumber
%
%\end{align}
%
%\end{small}
%where the $J_{i,j}(z)$ and $\bar{J}_{i,j}(z)$($i,j=1, 2,..., 5$) are the nonlinear term evaluated at steady-state by ignoring the 2nd order and higher order terms, such as $J_{1,1}=-A_1\rho_{e,ss}e^{-\frac{E_1}{T_{c,ss}}}$, $J_{1,2}=A_1(\rho_{s1}^{0}-\rho_{s1,ss})e^{-\frac{E_1}{T_{c,ss}}}$.
%The input operator $\mathcal{B}_1$ is defined as a bounded operator $\mathcal{B}_1$ =[ 0; 0; 0; 0; $b(z)$] and $b(z)=\frac{1}{2\bar{z} _{L}}1_{[z _{L}-\bar{z} _{L},z _{L}+\bar{z} _{L}]}(z)$. The controlled output operator $\mathcal{C}_2$ is determined as $\mathcal{C}_2(\cdot) =\text{diag}\{\int_{0}^{L}\delta (z -L) (\cdot )d\eta, \!0, \!0, \!0, \!0\}$, \! \!and the measured output operators are defined as $\mathcal{C}_{m1} \! \!= \! \!\text{diag}\{\int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta\}$ and $\mathcal{C}_{m2}\!=\!\text{diag}\{\int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z ) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z ) (\cdot )d\eta\}$, respectively.
%In this form, the domains of operator $\mathcal{A}_1$ and $\mathcal{A}_2$ are $\mathcal{D}(\mathcal{A}_1) \!= \! \{\left [ \phi _1 (z); \phi _2 (z); \phi _3 (z); \phi _4 (z); \phi _5 (z)\right ] \!\in L_{2}(0,1)^5|\phi _i (z)\text{ are abs. cont.,} \text{ and }  \phi _{i}(0)\!=\!0, i=1,2,...,5 \} $, and $\mathcal{D}(\mathcal{A}_2)\!=\! \{ [ \phi _6 (z); \phi _7 (z); \phi _8 (z); \phi _9 (z);\\ \phi _{10} (z) ]\in L_{2}(0,1)^5|\phi _j (z)\text{ are abs. cont.,}  j=6,7,...,10 \text{ and } \phi _{j}(L)=0, j=8,10\} $, respectively.
%Mathematically, the dynamics of the above cascade system (cook zone and wash zone) can be described by the following extended system:
%\begin{equation}
%\begin{aligned}
%\dot{x}(z,t)&=\mathcal{A}x(z,t)+\mathcal{B}u(t)\\
%y_c(t)&=\mathcal{C}_cx(z,t) \\
%y_m(t)&=\mathcal{C}_mx(z,t)
%\label{exsys} 
%
%\end{aligned}
%
%\end{equation}
%with boundary conditions given below:
%\begin{equation}
%\begin{aligned}
%x_{1} (0, t)&=0,
%x_{2} (0, t)=0,
%x_{3} (0, t)=0, 
%x_{4} (0, t)=0,
%x_{5} (0, t)=0\\
%x_{6} (0, t)&=x_{1} (L, t),
%x_{7} (0, t)=x_{2} (L, t),
%x_{8} (L, t)=0, 
%x_{9} (0, t)=x_{4} (L, t),
%x_{10} (L, t)=0
%\label{newbc} 
%
%\end{aligned}
%
%\end{equation}
%where the extended state $x(z,t)=[x^{(\mathcal{C})}(z,t);x^{(\mathcal{W})}(z,t)]$ and the corresponding extended operators are given as follows:
%\begin{equation}
%\begin{aligned}
%\mathcal{A}(\cdot)&=V\frac{\partial }{\partial z}(\cdot)+\psi(z)(\cdot), V=\text{bdiag}\{V_1,V_2\}, \psi (z)=\text{bdiag}\{\psi_1,\psi_2\}\\
%\mathcal{B}&=[\mathcal{B}_1;\textbf{0} ], \mathcal{C}_c=[\textbf{0} ,\mathcal{C}_2], 
%\mathcal{C}_m=\text{bdiag}\{\mathcal{C}_{m1},\mathcal{C}_{m2}\}
%
%\end{aligned}
%
%\end{equation}
%The state $x(z,t) \in \mathcal{X}$, with $\mathcal{X}={{L}}^2((0,L),{\mathbb{R}}^{10})$ being defined as a real separable Hilbert space with inner product $\langle\cdot, \cdot\rangle$. The input $u(t) \in {{L}}_{loc}^2([0,\infty),\mathcal{U})$ and output $y(t) \in {{L}}_{loc}^2([0,\infty),\mathcal{Y}) $, where $\mathcal{U}$ and $\mathcal{Y}$ are real separable Hilbert spaces.
%Based on the coupling conditions Equation (\ref{bcwash}), we have the domain of $\mathcal{A}$, which is $\mathcal{D}(\mathcal{A})=\mathcal{D}(\mathcal{A}_1)\bigoplus \mathcal{D}(\mathcal{A}_2) = \{[\phi _1;\phi _2;...;\phi _{10}]\in L_{2}(0,L)^{10}|[\phi _1;\phi _2;...;\phi _{5}] \in \mathcal{D}(\mathcal{A}_1), [\phi _6;\phi _7;...;\phi _{10}] \in \mathcal{D}(\mathcal{A}_2), \text{ and } \phi _{1}(L)=\phi _{6}(0),\phi _{2}(L)=\phi _{7}(0),\phi _{4}(L)=\phi _{9}(0) \} $ (missing citation). 
%The adjoint operator $\mathcal{A}^{*}$ is easily found using the inner product formula, $\langle\mathcal{A}\varphi_i, \phi_i \rangle = \langle\varphi_i, \mathcal{A}^*\phi_i \rangle, i=1,2,...,10$, and is:
%\begin{align}
%\mathcal{A}^{*}(\cdot)=-V\frac{\partial {(\cdot)}}{\partial z }+\psi^{*}(z)(\cdot)
%
%\end{align}
%with its domain defined as $ \mathcal{D}(\mathcal{A}^*)=\{[\phi _1;\phi _2;...;\phi _{10}]\in L_{2}(0,1)^{10},$
%$ \phi_i (z ) \text{ is absolutely continuous, } $\\$\frac{d \phi_i}{d z}\in L_{2}(0,1), \!\text{ with } i\!=\!1,2,...,10, \!\text{ and }\phi _{j}(L)\!=\!0, j=3,5,6,7,9,\phi _{m}(0)\!=\!0, m=\!8,10, \phi _{1}(L)\!=\!\frac{v_6}{v_1}\phi _{6}(0), \phi _{2}(L)=\frac{v_7}{v_2}\phi _{7}(0), \phi _{4}(L)=\frac{v_9}{v_4}\phi _{9}(0)\}$.
%\section*{Model Discretization}
%In this section, the Cayley-Tustin discretization framework is considered and applied to the linearized digester system without any spatial approximation induced. By the use of Cayley–Tustin transformation, a discrete-time state-space model for describing the extended digester system is established and realized by determining the resolvent operator, which is amenable to the discrete observer and controller designs. Meanwhile,  the essential properties of the continuous-time system stay invariant under this transformation, such as conservative (energy preserving)(missing citation), stability(missing citation), observability(missing citation); (missing citation), controllability.
%\subsection*{Cayley-Tustin Time Discretization Framework}
%The above linear system in Equation (\ref{exsys}) is considered. For a given time discretization $h>0$, and for $j\geq 1$ the Cayley-Tustin discretization is given by:
%\begin{equation}
%\begin{aligned}
%\frac{{x(jh) - x((j - 1)h)}}{h} &\approx \mathcal{A}\frac{{x(jh) + x((j - 1)h)}}{2} + \mathcal{B}u(jh)\\
%y_c(jh) &\approx \mathcal{C}_c\frac{{x(jh) + x((j - 1)h)}}{2} \\
%y_m(jh) &\approx \mathcal{C}_m\frac{{x(jh) + x((j - 1)h)}}{2} 
%
%\end{aligned}
%\label{cayley} 
%
%\end{equation}
%with $x(0) = {x_0}$, where the spatial dependence of $x$ is omitted for brevity. Then, let $\frac{u_j^{(h)}}{\sqrt{h}}$ be an approximation of $u(jh)$ by the mean value within a given sampling time, $\frac{u_j^{(h)}}{\sqrt{h}}=\frac{1}{h}\int_{(j-1)h}^{jh}u(t)dt$. It has been shown in (missing citation) that $\frac{u_j^{(h)}}{\sqrt{h}}$ converges to $u(jh)$ as $h \rightarrow 0$ in several different ways, similar for $Y(jh)$. Further, rewriting Equation (\ref{cayley}) gives the discrete time dynamics Equation (\ref{tustin}). It is frequently called Tustin discretization in the engineering literature, which is proposed in 1940s by Tustin and referred as Tustin transform in digital and sample-data control literature (missing citation).
%\begin{equation}
%\begin{aligned}
%\frac{{x_j- x_{j-1}}}{h} &{\approx \mathcal{A}\frac{{x_j+ x_{j-1}}}{2} + \mathcal{B}\frac{u_j}{\sqrt{h}}}\\
%{\frac{y_{cj}}{\sqrt{h}} }&{\approx \mathcal{C}_c\frac{{x_j + x_{j-1}}}{2} }\\
%{\frac{y_{mj}}{\sqrt{h}} }&{\approx \mathcal{C}_m\frac{{x_j + x_{j-1}}}{2} }
%
%\end{aligned}
%\label{tustin} 
%
%\end{equation}
%Through some basic computations, the following infinite-dimensional discrete-time state space model is obtained:
%\begin{equation}
%\begin{split}
%x_j &= {\mathcal{A}_d}x_{j-1}+ {\mathcal{B}_d}u_j \\
%y_{cj} &= {\mathcal{C}_{cd}}x_{j-1} + {\mathcal{D}_{cd}}u_j\\
%y_{mj}&= {\mathcal{C}_{md}}x_{j-1}+ {\mathcal{D}_{md}}u_j
%
%\end{split}
%
%\end{equation}
%where $\mathcal{A}_d$, $\mathcal{B}_d$, $\mathcal{C}_d$, $\mathcal{D}_d$, $ \mathcal{C}_{cd}$ and  $\mathcal{C}_{md}$ are the discrete-time spatial operators and we denote:
%\begin{align}
%\begin{pmatrix}
%\!\!\mathcal{A}_d &  \mathcal{B}_d\\ 
% \mathcal{C}_{cd}& \mathcal{D}_{cd} \\
%  \mathcal{C}_{md}& \mathcal{D}_{md}
%\end{pmatrix}
%\!\!=
%\begin{pmatrix}
%-I+2\delta \mathcal{R}(z,\delta ) & \sqrt {2\delta}\mathcal{R}(z,\delta ) \mathcal{B}\\ 
%\sqrt {2\delta}\mathcal{C}_c\mathcal{R}(z,\delta ) & \mathcal{G}_c(\delta )\\
%\sqrt {2\delta}\mathcal{C}_m\mathcal{R}(z,\delta ) & \mathcal{G}_m(\delta )
%\end{pmatrix}
%\label{operators} 
%
%\end{align}
%where $\delta=2/h$ and the resolvent is $\mathcal{R}(z,\delta)=(\delta -\mathcal{A})^{-1}$. Clearly, one must satisfy $\delta \in \rho(\mathcal{A})$ so that the resolvent operator is well-defined. In particular, $\mathcal{G}_c(\delta )$ denotes the transfer function from input to controlled output $\mathcal{G}_c(\delta )=\mathcal{C}_c(\delta  - \mathcal{A})^{ - 1}\mathcal{B}$, and $\mathcal{G}_m(\delta )$ denotes the transfer function from input to measured outputs $\mathcal{G}_m(\delta )=\mathcal{C}_m(\delta  - \mathcal{A})^{ - 1}\mathcal{B}$.
%The unbounded operators $\mathcal{A} $ of the continuous-time system are mapped into bounded operators $\mathcal{A}_d$ in the discrete-time counterpart through Cayley transform, which describe by the resolvent operator with $\mathcal{A}_d(\cdot )=[\delta  -\mathcal{A}]^{-1}[\delta  +\mathcal{A}](\cdot )=-I(\cdot )+2\delta[\delta  -\mathcal{A}]^{-1}(\cdot )=-I(\cdot )+2\delta\mathcal{R} (\delta,\mathcal{A})(\cdot )$. In addition, it has been demonstrated that the controllability and stability are invariant under this transformation (missing citation).
%\subsection*{Resolvent operator}
%The resolvent operator can be obtained by utilizing the Laplace transform. Under the zero input condition, we can have the following expression,
%\begin{align}
%\frac{\partial x(z,s)}{\partial z}=V^{-1}(sI-\psi)x(z,s)-V^{-1}x(z,0)
%
%\end{align}
%where $\psi$ is taken as the average of $\psi(z)$. By solving the above ODE, one obtains:
%\begin{align}
%x(z,s)=e^{\{V^{-1}(sI-\psi)z\}}x(0,s) -\int_{0}^{z}e^{\{V^{-1}(sI-\psi)(z-\eta)\}}V^{-1}x(\eta ,0)d\eta
%\label{lap} 
%
%\end{align}
%Since $\psi$ is the block diagonal matrix and $V$ is already in diagonal form, the matrix $F$ is finally a block diagonal matrix, denoted as follows:
%\begin{align}
%F=V^{-1}(sI-\psi)=
%\begin{bmatrix}
%[F_1]_{5\times 5} & \mathbf{0}_{5\times 5}\\ 
%\mathbf{0}_{5\times 5} & [F_2]_{5\times 5}
%
%\end{bmatrix}
%
%\end{align}
%Then $e^{Fz}=\text{diag}(e^{F_1z}, e^{F_2z})$ and we denote $[M_{ij}(z,s)]_{10\times 10}=e^{Fz }$ for simplicity.
%In order to obtain the unknown boundary conditions, one needs to evaluate the boundary conditions in Equation (\ref{newbc}) as follows:
%(1) At $z = 0$, one can plug $x_{1} (0, t)=0, x_{2} (0, t)=0, x_{3} (0, t)=0,  x_{4} (0, t)=0, x_{5} (0, t)=0$ into Equation (\ref{lap}) which leads to $M_{ij}(0,s)=0, i=1,2,...5, j=6,7,...10$.
%(2) At $z = L$, one can firstly substitute $x_{6} (0, s)=x_{1} (L, t), x_{7} (0, s)=x_{2} (L, t),x_{9} (0, s)=x_{4} (L, t)$ into Equation (\ref{lap}) which yields
%\begin{equation}
%\begin{aligned}
%x_{6} (0, s)=-\sum_{k=1}^{10}\int_{0}^{L}M_{1k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \\
%x_{7} (0, s)=-\sum_{k=1}^{10}\int_{0}^{L}M_{2k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \\
%x_{9} (0, s)=-\sum_{k=1}^{10}\int_{0}^{L}M_{4k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta
%
%\end{aligned}
%
%\end{equation}
%where $V_v$ denotes $V^{-1}$.
%Then, substituting the boundary conditions $x_{8} (L, s)=0,x_{10} (L, s)=0$, one can have
%\begin{small}
%\begin{align}
%&
%\begin{bmatrix}
% M_{88}&  M_{810}\\ 
%  M_{108}&  M_{1010}
%\end{bmatrix}
%\begin{bmatrix}
%x_8(0,s)\\ 
%x_{10}(0,s)
%\end{bmatrix}
%\nonumber\\ 
%&=
%\begin{bmatrix}
%\sum_{k=1}^{10}\int_{0}^{L}M_{8k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta-M_{86}x_6(0,s)-M_{87}x_7(0,s)-M_{89}x_9(0,s)\\ 
%\sum_{k=1}^{10}\int_{0}^{L}M_{10k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta-M_{106}x_6(0,s)-M_{107}x_7(0,s)-M_{109}x_9(0,s)
%
%\end{bmatrix}
%
%\end{align}
%
%\end{small}
%By substituting $x_{6} (0, s)$, $x_{7} (0, s)$, $x_{9} (0, s)$ and solving the above equations, we can obtain the expressions of $x_{8} (0, s)$ and $x_{10} (0, s)$ as follows,
%\begin{align}
%x_8(0,s)=&G_{11}(\sum_{k=1}^{10}\int_{0}^{L}M_{8k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta)+G_{12}(\sum_{k=1}^{10}\int_{0}^{L}M_{10k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta) \nonumber\\ 
%&+(G_{11}M_{86}+G_{12}M_{106})\sum_{k=1}^{10}\int_{0}^{L}M_{1k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \nonumber \\
%&+(G_{11}M_{87}+G_{12}M_{107})\sum_{k=1}^{10}\int_{0}^{L}M_{2k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \nonumber\\
%&+(G_{11}M_{89}+G_{12}M_{109})\sum_{k=1}^{10}\int_{0}^{L}M_{4k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \\
%x_{10}(0,s)=&G_{21}(\sum_{k=1}^{10}\int_{0}^{L}M_{8k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta)+G_{22}(\sum_{k=1}^{10}\int_{0}^{L}M_{10k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta) \nonumber\\ 
%&+(G_{21}M_{86}+G_{22}M_{106})\sum_{k=1}^{10}\int_{0}^{L}M_{1k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \nonumber\\
%&+(G_{21}M_{87}+G_{22}M_{107})\sum_{k=1}^{10}\int_{0}^{L}M_{2k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \nonumber\\
%&+(G_{21}M_{89}+G_{22}M_{109})\sum_{k=1}^{10}\int_{0}^{L}M_{4k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta
%
%\end{align}
%where
%\begin{align}
%G=
%\begin{bmatrix}
% M_{88}&  M_{810}\\ 
%  M_{108}&  M_{1010}
%\end{bmatrix}
%^{-1}
%
%\end{align}
%Therefore, the resolvent operator is determined as following form by substituting the boundary conditions (Equation (\ref{newbc})):
%\begin{align}
%x(z,s)=[\mathcal{R}_{ij}(z,s)]_{10\times 10}x(\eta,0)
%
%\end{align}
%where
%\begin{align}
%\mathcal{R}_{ij}=\left\{
%\begin{matrix}
%\mathcal{R}c_{ij}, i\in [1,5], j\in [1,10]\\ 
%\mathcal{R}w_{ij}, i\in [6,10], j\in [1,10]
%\end{matrix}
%\right.
%
%\end{align}
%where
%\begin{align}
%\mathcal{R}c_{ij}=-\int_{0}^{z}M_{ij}(z-\eta ,s)V_{vjj}(\cdot )d\eta, i\in [1,5], j\in [1,10] 
%
%\end{align}
%\begin{small}
%\begin{equation}
%\begin{aligned}
%\mathcal{R}w_{ij}&=(M_{i8}(z,s)G_{11}+M_{i10}(z,s)G_{21})\int_{0}^{L}M_{8j}(L-\eta,s )V_{vjj}(\cdot )d\eta \nonumber\\
%&+(M_{i8}(z,s)G_{12}+M_{i10}(z,s)G_{22})\int_{0}^{L}M_{10j}(L-\eta,s )V_{vjj}(\cdot )d\eta \nonumber\\
%&+(M_{i8}(z,s)G_{11}M_{86}(L,s)+M_{i8}(z,s)G_{12}M_{106}(L,s)+M_{i10}(z,s)G_{21}M_{86}(L,s) \nonumber\\
%&+M_{i10}(z,s)G_{22}M_{106}(L,s)-M_{i6}(z,s))\int_{0}^{L}M_{1j}(L-\eta,s )V_{vjj}(\cdot )d\eta
%  \nonumber\\
%&+(M_{i8}(z,s)G_{11}M_{87}(L,s)\!+\!M_{i8}(z,s)G_{12}M_{107}(L,s)\!+\!M_{i10}(z,s)G_{21}M_{87}(L,s)\!+\!M_{i10}(z,s)G_{22}M_{107}(L,s) \nonumber\\
%&-M_{i7}(z,s))\int_{0}^{L}M_{2j}(L-\eta,s )V_{vjj}(\cdot )d\eta+(M_{i8}(z,s)G_{11}M_{89}(L,s)+M_{i8}(z,s)G_{12}M_{109}(L,s) \nonumber\\
%&+M_{i10}(z,s)G_{21}M_{89}(L,s)+M_{i10}(z,s)G_{22}M_{109}(L,s)-M_{i9}(z,s))\int_{0}^{L}M_{4j}(L-\eta,s )V_{vjj}(\cdot )d\eta \nonumber\\
%&-\int_{0}^{z}M_{ij}(z-\eta ,s)V_{vjj}(\cdot )d\eta, i\in [6,10], j\in [1,10]
%
%\end{aligned}
%
%\end{equation}
%
%\end{small}
%The discrete-time operators in Equation (\ref{operators}) can be solved by straightforwardly substituting the above resolvent operators. Afterwards, the discrete-time linear model is obtained:
%\begin{equation}
%\begin{split}
%x(z ,k) &= {\mathcal{A}_d}x(z ,k - 1) + {\mathcal{B}_d}u(k)\\
%y_c(k) &= {\mathcal{C}_{cd}}x(z ,k - 1) + {\mathcal{D}_{cd}}u(k)\\
%y_m(k) &= {\mathcal{C}_{md}}x(z ,k - 1) + {\mathcal{D}_{md}}u(k)
%
%\end{split}
%\label{discretetimemodel} 
%
%\end{equation}
%with the boundary conditions (Equation (\ref{newbc})).
%\section*{Observer-Based MPC Design}
%An observer-based model predictive controller is designed for the discrete-time pulp digester system. In particular, a discrete Luenberger observer is designed first to reconstruct the states based on the available real-time measurements. The Luenberger observer is one of the practical and easy-to-realize observer, which is further considered in a discrete setting controller realization.
% The constrained optimal controller design for the finite-dimensional system theory is extended and deployed for the infinite-dimensional digester system. The overall closed-loop operation of the digester process is schematically presented in Figure \ref{framework}.
%\subsection*{Discrete Luenberger observer design}
%Firstly, let us recall the linearized discrete-time model:
%\begin{equation}
%\begin{split}
%x(z ,k) &= {\mathcal{A}_d}x(z ,k - 1) + {\mathcal{B}_d}u(k)\\
%y_m(k) &= {\mathcal{C}_{md}}x(z ,k - 1) + {\mathcal{D}_{md}}u(k)
%
%\end{split}
%\label{dtmodel} 
%
%\end{equation}
%The discrete Luenberger observer is presented by the following standard form:
%\begin{equation}
%\begin{split}
%\hat{x}(\zeta ,k) &= {\mathcal{A}_d}\hat{x}(\zeta ,k - 1) + {\mathcal{B}_d}u(k) +L_d(y_m(k)-\hat{y}_m(k)) \\
%\hat{y}_m(k)&={\mathcal{C}_{md}}\hat{x}(\zeta ,k - 1) + {\mathcal{D}_{md}}u(k)
%
%\end{split}
%\label{disluen} 
%
%\end{equation}
%where the reconstructed state $\hat{x}(\zeta ,k)$ is defined as a copy of
%the system dynamics and $L_d$ is the discrete observer gain to be designed. Stability of the observer implies that the state estimation error, $e_k=x(\zeta, k)-\hat{x}(\zeta, k)$, converges to zero within a certain time. The error dynamic equation is shown as follows:
%\begin{align}
%e_k=(\mathcal{A}_d-L_d\mathcal{C}_d)e_{k-1}=\tilde{\mathcal{A}}_de_{k-1}
%\label{errordy} 
%
%\end{align}
%To guarantee the operator $\tilde{\mathcal{A}}_d$ in the state
%estimation error dynamics given by Equation (\ref{errordy}) is stable, the design objective is to determine the appropriate spatially varying gain $L_d$.
%By Curtain and Zwart (missing citation), it can be shown that the operator $\tilde{\mathcal{A}}_d$ is power stable if and only if there exists a non-negative self-adjoint operator $Q_d$ such that
%\begin{align}
%\tilde{\mathcal{A}}_dQ_d\tilde{\mathcal{A}}_d^*-Q_d=- M_d, \text{on }    \mathcal{X}
%
%\end{align}
%where $M_d$ is a positive definite design parameter.
%Then, let us assume that the pair $(\mathcal{A}, \mathcal{C})$ is exponential detectable, then if there exists a nonnegative self-adjoint operator $Q_d$ which is the solution of the following operator Riccati equation (missing citation); (missing citation):
%\begin{align}
%\mathcal{A}_dQ_d\mathcal{A}_d^*-Q_d-L_d(P+2I)L_d^*+M_d=0, \text{on }    \mathcal{X}
%\label{disricca} 
%
%\end{align}
%where
%\begin{align}
%L_d=T(P+I)^{-1},
%T=\mathcal{A}_dQ_d\mathcal{C}_d^*,
%P=\mathcal{C}_dQ_d\mathcal{C}_d^*  \nonumber
%
%\end{align}
%The observer gain $L_d= \mathcal{A}_dQ_d\mathcal{C}_d^*(\mathcal{C}_dQ\mathcal{C}_d^*+I)^{-1}$ is a strongly stabilizing gain which guarantees the power stability of $\tilde{\mathcal{A}}_d=\mathcal{A}_d-L_d\mathcal{C}_d$. To solve the algebraic Riccati Equation (\ref{disricca}), one can utilize the numerical iteration methods, such as Newton-Kleinman iteration method (missing citation), and the detailed procedures were provided (missing citation).
%\subsection*{MPC design for target tracking}
%The ultimate objective of a pulp and paper mill is to ensure the specified quality of the end products while meeting the production targets and minimizing the operational costs. Specifically, when it comes to practical applications, the most characteristic point is that the species of feed wood chips are switched frequently, softwood and hardwood chips, to supply the required amount of pulp according to the production schedule of the paper machines (missing citation). Accordingly, the desired product quality needs to change in the middle of the operation, i.e., the shift operations. To realize it, the MPC is developed for the infinite-dimensional setting, emerging from the finite-dimensional linear time-invariant systems,
%see Rawlings et al. (missing citation).
%In this case, we consider that the system output is required to track a nonzero target vector, $y_t$, then state and input vectors, $x_t$ and $u_t$, are required which bring the system to $y_t$ at steady-state (missing citation); (missing citation); (missing citation). The state and input target can be computed by solving the following quadratic program.
%\begin{equation}
%\begin{split}
%\mathop {\min }\limits_{{x_t, u_t}} \;  & { (u_t - \bar{u}) ^TR_t(u_t - \bar{u}) }\\
%{\rm{s.t. \;\;}}&
%\begin{bmatrix}
%I-A_d & -B_d\\ 
%C_d & D_d
%\end{bmatrix}
%\begin{bmatrix}
%x_t\\ 
%u_t
%\end{bmatrix}
%=
%\begin{bmatrix}
%0\\ 
%y_t
%\end{bmatrix}
% \\
%&{u_{\min }} \le u_t \le {u_{\max }} \\
%&{y_{\min }} \le y_t \le {y_{\max }} 
%
%\end{split}
%\label{stateinput} 
%
%\end{equation}
%In this quadratic program, $\bar{u}$ is the set point for the manipulated variables and $R_t$, is symmetric positive definite. Notice that, often, the input set point is not specified and it can be assumed zero in order to use (missing citation). The equality constraints guarantee a steady-state solution and offset free tracking of the target vector (missing citation). $u_{min}$, $u_{max}$, $y_{min}$, $y_{max}$ are the input and output constraints, respectively.
%Then, the the following quadratic objective function is used for the regulator to track a nonzero target vector.
%\begin{equation}
%\begin{split}
%\mathop {\min }\limits_{{u^N}} \;  &\sum\limits_{j = 0}^{N - 1}\;  { {\langle ( y_{k + j|k}-y_t ) ,Q(y_{k + j|k} -y_t) \rangle}  + \langle (u_{k + j|k}-u_t), R(u_{k + j|k}-u_t)\rangle} \\
%& + \langle (x_{k + N-1|k}-x_t), \bar Q(x_{k + N-1|k}-x_t)   \rangle  \\
%{\rm{s.t. \;\;}}&\text{Lunberger observer, Equation }(\ref{disluen})\\
%&{u_{\min }} \le u_{k + j|k} \le {u_{\max }} \\
%&{y_{\min }} \le y_{k + j|k} \le {y_{\max }} 
%
%\end{split}
%\label{trackmpc} 
%
%\end{equation}
%where $N$ is the prediction horizon, and $Q$, $R$ are symmetric positive semidefinite and symmetric positive definite spatial operator, respectively. $ \bar{Q}$ is the spatial operator to penalize the terminal state which depends on the stability of the given model. 
%The target vector $x_t$, and $u_t$ are computed from the quadratic program in Equation (\ref{stateinput}). $y_{k + j|k}$ and $u_{k + j|k}$ represent the output and input variables at future time $k+j$ predicted at current time $k$, and the term $\Delta u_{k + j|k}$ denotes the change of an input vector at time $k+j$ as $\Delta u_{k + j|k} = u_{k + j|k}-u_{k + j-1|k}$. The vector $u^N$ includes the control sequence $\{u_{k|k}, u_{k + 1|k}, u_{k + 2|k}, ..., u_{k + N-1|k}\}$ and the first element $u_{k|k}$ will be injected to the plant as the future control action.
%As discussed by H. Kwakernaak and R. Sivan (missing citation), using the targets computed from Equation (\ref{stateinput}), we define a shifted input $\tilde{u}_k=u_k - u_t$, a shifted state $\tilde{x}_k=\hat{x}_k - x_t$, and a shifted output $\tilde{y}_k=\hat{y}_k - y_t$ to reduce the problem to the standard form. The corresponding constraints can be translated to constraints on $\tilde{u}$ and $\tilde{y}$.
%Thus the regulator optimization problem Equation (\ref{trackmpc}) becomes:
%\begin{equation}
%\begin{split}
%\mathop {\min }\limits_{{u^N}} \;  &\sum\limits_{j = 0}^{N - 1}\;  { {\langle \tilde{y}_{k + j|k} ,Q\tilde{y}_{k + j|k} \rangle}  + \langle \tilde{u}_{k + j|k}, R\tilde{u}_{k + j|k}\rangle}  + \langle \tilde{x}_{k + j|k}, \bar Q\tilde{x}_{k + j|k}   \rangle  \\
%{\rm{s.t. \;\;}}&\tilde{x}_{k + j|k} = {\mathcal{A}_d}\tilde{x}_{k + j-1|k} + {\mathcal{B}_d}\tilde{u}_{k + j|k}  \\
%&\tilde{y}_{k + j|k} = {\mathcal{C}_d}\tilde{x}_{k + j-1|k} + {\mathcal{D}_d}\tilde{u}_{k + j|k} \\
%&{u_{\min }-u_{t|k}} \le \tilde{u}_{k + j|k} \le {u_{\max }-u_{t|k}} \\
%&{y_{\min }-y_{t|k}} \le \tilde{y}_{k + j|k} \le {y_{\max }-y_{t|k}} 
%
%\end{split}
%\label{trackmpc2} 
%
%\end{equation}
%According to the nature of transport reaction systems, the operator $ \bar{Q}$ can be determined from the positive self-adjoint solution of the following discrete-time Lyapunov equation:
%\begin{align}
%{\mathcal{A}^*_d}\bar Q \mathcal{A}_d  - \bar Q
%=- {\mathcal{C}^*_d}Q\mathcal{C}_d
%\label{ctlya} 
%
%\end{align}
%or equivalently the continuous-time Lyapunov equation (missing citation):
%\begin{align}
%{\mathcal{A}^*}\bar Q + \bar Q\mathcal{A} =  - {\mathcal{C}^*}Q\mathcal{C}
%\label{dtlya} 
%
%\end{align}
%on the dual space of $\mathcal{X}_{-1}$. In addition, the operator $ \bar{Q}$ is the unique positive self-adjoint solution of the Lyapunov equations (Equation (\ref{ctlya}) and (\ref{dtlya})) (missing citation); (missing citation).
%Before further manipulate the objective function (Equation (\ref{trackmpc2})), we introduce the following notations: $Y_k=\{y_{k+n}\}_{n=1}^{N}\in \mathcal{Y}^{N}$ and $U_k=\{u_{k+n}\}_{n=1}^{N}\in \mathcal{U}^{N}$. As a result, the straightforward algebraic manipulation of the objective function presented in Equation (\ref{trackmpc2}) leads to the following quadratic programming optimization problem:
%\begin{equation}
%\begin{split}
%\mathop {\min }\limits_{U_k} \;& \langle U_k, HU_k \rangle + 2\langle U_k,F\tilde{x}_k  \rangle  +   \langle \tilde{x}_k,\bar Q\tilde{x}_k  \rangle  \\
%{\rm{s.t. \;\;}}&{U_{\min }} \le U_k \le {U_{\max }} \\
%&{Y_{\min }} \le GU_k + S\tilde{x}_k \le {Y_{\max }}
%
%\end{split}
%\label{qprog} 
%
%\end{equation}
%where $H \in \mathcal{L}( \mathcal{U}^N)$ is positive and self-adjoint, which is given by:
%\begin{align}
%h_{i,j}=\left\{
%\begin{matrix}
%\mathcal{D}_{d}^{*}Q\mathcal{D}_d+\mathcal{B}_{d}^{*}\bar{Q}\mathcal{B}_d+R & \text{for }  i=j \\
%\mathcal{D}_{d}^{*}Q\mathcal{C}_d\mathcal{A}_{d}^{i-j-1}\mathcal{B}_d+\mathcal{B}_{d}^{*}\bar{Q}\mathcal{A}_{d}^{i-j}\mathcal{B}_d &  \text{for } i>j \\
%h_{j,i}^{*} &  \text{for } i<j
%\end{matrix}
%\right.
%
%\end{align}
%and $F$ is given by $F=\{\mathcal{D}_d^*Q{\mathcal{C}_d}\mathcal{A}_{d}^{k-1}+\mathcal{B}_d^*\bar Q\mathcal{A}_d^k\}_{k=1}^{N-1}$. The mateix $G$ is a lower triangular given by
%\begin{align}
%g_{i,j}=\left\{
%\begin{matrix}
%\mathcal{D}_{d} & \text{for }  i=j \\
%\mathcal{C}_d\mathcal{A}_{d}^{i-j-1}\mathcal{B}_d &  \text{for } i>j \\
%0 &  \text{for } i<j
%\end{matrix}
%\right. \nonumber
%
%\end{align}
%and $S=\{{\mathcal{C}_d}\mathcal{A}_{d}^{k-1}\}_{k=1}^N$.
%The inner products in the objective function given in Equation (\ref{qprog}) are vector products as $U$ is the finite-dimensional input space, and therefore we have a finite dimensional quadratic optimization problem:
%\begin{equation}
%\begin{split}
%\mathop {\min }\limits_{U_k} \;&U_k^THU_k  + 2U_k^TF\tilde{x}_k    \\
%{\rm{s.t. \;\;}}&
%\begin{bmatrix}
%I\\ 
%-I\\ 
%G\\ 
%-G
%\end{bmatrix}
%U_k \le
%\begin{bmatrix}
%U_{\max }\\
%-U_{\min }\\
%Y_{\max }-S\tilde{x}_k\\
%-Y_{\min }+S\tilde{x}_k
%
%\end{bmatrix}
%
%\end{split}
%\label{qprognew} 
%
%\end{equation}
%Here we neglect the term $\langle \tilde{x}_k,\bar Q\tilde{x}_k  \rangle$ as $\tilde{x}_k$ is the initial condition for step $k+1$ and cannot be affected by the control input. Thus, the optimal input trajectory ($U_k$) can be obtained as the solution of the feasible quadratic optimization problem (Equation (\ref{qprognew})) converges to zero.
%\section*{Numerical Simulations}
%In this section, the closed-loop performance of the proposed MPC framework is demonstrated. The temperature of free liquor flowing into the cook zone is selected as the manipulated input variable and the concentration of lignin is selected as the controlled output variable. In this case, we consider the target tracking of the output by using the proposed MPC. The resulting constrained optimization problems become quadratic programming problems which are solved using the MATLAB subroutine QuadProg. Both the control and prediction horizons are chosen to be 50 sampling periods. The sampling time is set to be 10 min and the internal spatial discretization is taken as 0.05 m.
%The values of all system parameters taken in the simulations are listed in Table \ref{tab1}. For the initial conditions of the dynamic system, we consider $x_{1}(z, 0)=0.067z, x_{2}(z, 0)=0.1646z, x_{3}(z, 0)=0.2377z, x_{4}(z, 0)=1.7073z, x_{5}(z, 0)=0.8661z, x_{6}(z, 0)=0.0336+2.87\times10^{-4}z, x_{7}(z, 0)=0.8232+0.0095z, x_{8}(z, 0)=0.0012z, x_{9}(z, 0)=8.5363+0.1767z$, and  $x_{10}(z, 0)=-0.086z$. In addition, $\bar{\zeta} _{L}=0.15$ is chosen for the input operator. 
% The selected parameter values for MPC implementation are listed in Table \ref{tab2} .
% 
%As illustrated in Figure \ref{x}, the open-loop states converge to their corresponding steady-states rapidly which indicates the original plant is intrinsically stable. The figures in the left column represent the evolution of the states in the cook zone, and the right column denotes the evolution of the states in the wash zone.
%Figure \ref{x}(a) and \ref{x}(f) shows how the concentration of lignin decreases smoothly down toward the end of the cooking zone where the reactions
%are stopped (or quenched) by displacement of the hot liquor with dilute wash liquor from below. Hence, no significant decrease occurs in the wash zone. 
%The concentration of alkali of entrapped liquor phase in the feed flow to the cook zone is about $23$ g/l (as shown in Figure \ref{x}(b)) and then is consumed giving a decreasing profile down toward the extraction. At the bottom of the wash zone, the concentration of alkali of entrapped liquor phase is about $8$ g/l (as shown in Figure \ref{x}(g)). A similar trend occurs for the alkali of the free liquor phase, as shown in Figures \ref{x}(c) and \ref{x}(h).
%Figure \ref{x}(d) and \ref{x}(e) shows the temperature profiles for the two phases (wood chip phase and free liquor phase) at cook zone, and the temperature profiles of them at wash zone shown in Figure \ref{x}(i) and \ref{x}(j). The temperature of wood chips rises due to the exothermic reactions and is also affected by the high temperature of the circulation liquor. Below the extraction screens, in the wash zone, the chips are rapidly cooled down by the wash water which has a temperature of $360$K at the inlet in the bottom.
%By implementing the proposed MPC frameworks, we aim to steer the system to the desired targets without violating the physical constraints of actuators and sensors. In this case, we consider the shifted output targets, which are chosen as $y_{t1}=0.06$ and $y_{t2}=0.14$, respectively.
%First of all, the steady-state target $x_t$ and input vector $u_t$ are computed from the quadratic program Equation (\ref{stateinput}). Then, based on the operator Riccati Equation (\ref{disricca}), we determine $Q_d$ in the discrete-time observer design. The estimate state $\hat{x}$ from the Luenberger observer is then utilized for the MPC design. Finally, the optimal input trajectory is computed by solving the above optimization problem outlined in Equation (\ref{qprognew}) with a receding prediction horizon strategy.
% The state evolutions of the closed-loop system under the MPC law are obtained and shown in Figure \ref{closex}. Comparing with the open-loop state profiles, the closed-loop system is able to track the target steady-states when the grade shift happens ($t=370$ min).
%In addition, three pairs of states in cook zone and wash zone are successfully connected through the extend system (Equation (\ref{exsys})), such as, $\rho_{s}^{\mathcal{C}}(z,t)$ and $\rho_{s}^{\mathcal{W}}(z,t)$; $\rho_{e}^{\mathcal{C}}(z,t)$ and $\rho_{e}^{\mathcal{W}}(z,t)$; $T_c^{\mathcal{C}}(z,t)$ and $T_c^{\mathcal{W}}(z,t)$.
%The open-loop output profile and the close-loop profile under the observer-based MPC law are shown in Figure \ref{tr}. Without implementing the controller, the concentration of lignin converges to its steady-state, while violating the given constraints of the system. In the closed-loop system, the concentration of lignin is able to track the targets or desired values and satisfies the requirements of the constraint simultaneously. As the targets are switched, the output can also achieve target tracking across the original steady-state.
%Specifically, the target tracking above the steady-state is realized in the first period, $t\in[0,370]$, where $y_{t1}=0.06$ is considered. Similarly, the target tracking blows the steady-state is realized in the second period (i.t.,$t\in(370,740]$) when $y_{t2}=0.14$ is taken into account.
%Typically, this results in higher yield, thus lowering the operating cost significantly (missing citation).
%The free-liquor temperature profile, computed by the proposed model-based MPC system at each sampling time, is presented in Figure \ref{u}. 
%The input trajectory corresponds to the output variables, that is, in the first 370 minutes, the optimal input variables fluctuate between 425-430 and converge to 427K when the output goes to track target 1. In the last 370 minutes, the corresponding input rises and stabilizes at 436K when target 1 switches to target 2. In these two stages of tracking, the input variables are constrained within the given bounds of actuators. The simulation studies demonstrate that the extended system is able to describe the dynamics of the original cascade system which contains a co-current zone and a counter-current zone. Furthermore, it is possible to realize the optimal control of the output in the wash zone through the operation of the input in the cook zone. The effectiveness can be demonstrated from the proposed MPC design.
%\section*{Conclusions}
%In this work, dynamic modelling and model predictive control design of a continuous pulp digester described by ten linearized first-order coupled hyperbolic equations was developed. The connected cook zone and wash zone of the digester were modelled as a cascade PDE system. 
%By using Cayley–Tustin transformation, the linearized continuous-time infinite-dimensional model was transformed into a discrete-time infinite-dimensional model without spatial discretization and model reduction which preserves the input-output stability of the system. A Luenberger observer was designed to realize the state estimation of the system and the discrete-time Riccati equation was used to calculate the observer gain. The model predictive controller was formulated on that basis to realize target tracking and account for input and output constraints when it comes to the shift operations of the digester. The closed-loop simulation results have demonstrated that the controlled variables were able to reach to the target values and satisfy the constraints of actuators simultaneously.
%\newpage
%\newpage
%\newpage
%\listoffigures
%\newpage
%\newpage
%\newpage
%\listoftables
%\newpage
%\newpage
%\vspace{-1.6mm}
%\citestyle{C}
%\newtheorem{remark}{Remark}
%\newtheorem{lemma}{Lemma}
%\newtheorem{proposition}{Proposition}
%\newtheorem{corollary}{Corollary}
%\newtheorem{theorem}{Theorem}
%\newtheorem{alg}{Algorithm}
%\allowdisplaybreaks


\begin{document}

\title{Dynamic Modelling and Model Predictive Control of a Continuous Pulp
Digester}



\author[1]{Lu Zhang}%
\author[2]{Junyao Xie}%
\author[2]{Stevan Dubljevic}%
\affil[1]{Affiliation not available}%
\affil[2]{University of Alberta}%


\vspace{-1em}



  
  \date{May 24, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
%<<<<<<< HEAD
%This work explores the model predictive controller design of the continuous pulp digester process consisting of co-current zone and counter-current zone that are modeled by a set of nonlinear coupled hyperbolic partial differential equations (PDE). The system of interest is a not spectral and therefore all the arguments of low-dimensional representation by having just a few eigenmodes too describe dynamics do not hold and cannot be applied to in the given setting. To address this challenge, the nonlinear continuous-time model is linearized by applying the Taylor series expansion and then discretized utilizing the Cayley-Tustin discretization framework.
%This work explores the model predictive controller design of continuous pulp digester process consisting of co-current zone and counter-current zone modeled by a set of nonlinear coupled hyperbolic partial differential equations (PDE). The distributed parameter system of interest is a not spectral and low-dimensional representation designs by accounting for slow dynamics to describe model evolution does not hold and cannot be applied in the given setting. To address this challenge, the nonlinear continuous-time model is linearized and time discretized utilizing the Cayley-Tustin discretization framework.
%This approach guarantees the system theoretic properties and ensures energy and structure preservation without spatial discretization or model reduction. The discrete PDE model is utilized in the Kalman filter designed in order to estimate the spatial states and reject the noise. Further, a discrete model is used in model predictive controller design to regulate the output. Finally, the results of the state estimation are discussed, and the performance of controller is demonstrated by simulation studies.

%Specifically, in order to reduce the model complexity, we propose a simplified PDE model to describe the dominant dynamics of the delignification process. 
%A linearized model is obtained by applying the Taylor series expansion, which is then discretized utilizing the Cayley-Tustin discretization framework.
% In particular, this approach guarantees the system theoretic properties and ensures energy and structure preservation without spatial discretization or model reduction.  A Kalman filter is developed to estimate the states and reject the noise. Based on this, a model predictive controller is designed to regulate the output. Finally, the results of the state estimation are discussed, and the performance of controller is demonstrated by simulation studies.

%=======
This work explores the model predictive controller design of the continuous pulp digester process consisting of the co-current zone and counter-current zone modelled by a set of nonlinear coupled hyperbolic partial differential equations (PDE). 
%The distributed parameter system of interest is a not spectral and low-dimensional representation design by accounting for slow dynamics to describe model evolution does not hold and cannot be applied in the given setting. 
The distributed parameter system of interest is not spectral and slow-fast dynamic separation does not hold. 
To address this challenge, the nonlinear continuous-time model is linearized and discretized in time utilizing the Cayley-Tustin discretization framework, which ensures system theoretic properties and structure preservation without spatial discretization or model reduction. 
%Specifically, in order to reduce the model complexity, we propose a simplified PDE model to describe the dominant dynamics of the delignification process. 
%A linearized model is obtained by applying the Taylor series expansion, which is then discretized utilizing the Cayley-Tustin discretization framework. 
%This approach guarantees the system theoretic properties and ensures energy and structure preservation without spatial discretization or model reduction. 
%The discrete PDE model is utilized in the Luenberger observer designed in order to estimate the spatial states. Further, in order to realize shift operations of the digester, a discrete-time model is used in model predictive controller design to regulate the output and respect the system constraints. 
%The discrete model is used in full state model predictive controller design, which is augmented by the Luenberger observer design to realize shift operations of output and respect the system constraints. 
The discrete model is used in the full state model predictive controller design, which is augmented by the Luenberger observer design to achieve the output constrained regulation.
%The discrete PDE model is utilized in the Kalman filter designed in order to estimate the spatial states and reject the noise. Further, a discrete model is used in model predictive controller design to regulate the output. 
Finally, a numerical example is provided to demonstrate the feasibility and applicability of the proposed controller designs.
%Finally, the performance of controller and observer are demonstrated by simulation studies.
%>>>>>>> 6247784326c6d28dc13320b3d1779e933ecb2eb4


%, and the Taylor series expansion is applied for the model linearization. 
%
%In addition, the comparison between the softwood chip and hardwood chip in the digester are considered.
%
%The objective of this work is to realize the Kappa number tracking when the measurements are not directly available.%
\end{abstract}%



\sloppy


\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/paper/paper}
\caption{{\label{paper} The categories of paper products%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/digester/digester}
\caption{{\label{digester} Simplified scheme of a continuous pulp digester \protect(missing citation)%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/newsche/newsche}
\caption{{\label{scheme} The cook zone and wash zone in a digester%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/phaserepresent/phaserepresent}
\caption{{\label{phase} The conceptual model of the mass in a digester%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/ssopen/ssopen}
\caption{{\label{ss} Steady-state profiles of the digester (The solid lines denote the steady-states of the components in cook zone; the dash lines denote the steady-states of the corresponding components in wash zone.)%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/framenew/framenew}
\caption{{\label{framework} The proposed closed-loop operation framework%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/openx2/openx2}
\caption{{\label{x} The open-loop state profiles of the digester%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/closex5/closex5}
\caption{{\label{closex} The state profiles of the digester under closed-loop operation%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/trackobsvstair2/trackobsvstair2}
\caption{{\label{tr} The concentration profile of lignin under closed-loop operation%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/tracku2/tracku2}
\caption{{\label{u} The optimal manipulated input trajectory under closed-loop operation%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{table}[H]
\caption {{Notations and values of parameters }} \label{tab1} 
\vspace{1.6mm}
\centering
%\begin{small}
\begin{tabular}{c c c c c c c c c }
\hline
 \hline
Process parameters& Notations & Numerical Values \\
\hline
Volumetric flow rate of chip     & $\dot{V}_c $        & 0.0267 $m^3/min      $  \\
%(1.3964)             \\ 
Volumetric flow rate of free liquor      & $\dot{V}_f  $    & 0.09 $ m^3/min         $\\              
%(2.3497)\\ 
Digester cross sectional area     &$ A    $     & 1 $ m^2      $                     \\
Concentration of non-reactive lignin &$ \rho_{s}^{0}    $     & 0.015 $ kg/m^3      $                     \\
Frequency factor of lignin reactions &$ A_1    $     & 0.09 $ m^3/kg\cdot min      $                     \\ 
Activation energy for lignin &$ E_1    $     & 38 $ kJ/mol \cdot K      $                     \\ 
Stoichiometric coefficient for lignin reactions &$  b_1   $     & 0.15                     \\ 
Stoichiometric coefficient for carbohydrates reactions &$  b_2   $     & 0.25                     \\ 
Heat capacities of the wood &$C_{ps}    $     & 1.47 $  kJ/kg\cdot K      $                     \\
Heat capacities of the liquor &$C_{pl}    $     & 4.19 $  kJ/kg\cdot K      $                     \\
Heat of reaction&$\Delta H_R    $     & -581 $  kJ/kg      $                     \\
Water density&$\rho_w    $     & 1000 $  kg/m^3      $                     \\
 \hline
  \hline
\end{tabular}
%\end{small}

\end{table}\selectlanguage{english}
\begin{table}[H]
\caption {{Parameters for the MPC Design }} \label{tab2} 
\vspace{1.6mm}
\centering
\begin{tabular}{c c c}
\hline
 \hline
Descriptions       & Notations                   & Values            \\
\hline
Sampling time      & $h $                          & 10 $min$            \\
Prediction horizon & $N $                          & 50                \\
Input weight       & $R$                           & 0.5               \\
Output weight      & $Q$                           & 0.5               \\
Input constraints   & $[u_{min},u_{max}]$ & $[414.5K, 440.5K]$  \\
Controlled output constraints  &$ [y_{min},y_{max}{]}$ & $[0.04\%, 0.26\%]$ \\
 \hline
  \hline
\end{tabular}

\end{table}

\setstretch{1.5}

\newpage

\section*{Introduction}
The pulp and paper industry has a profound influence on the economy of the world, which produces pulp, paper, paperboard, and various cellulose-based products.
Even though electronic media and paperless communication have been widely expanded, the global pulp and paper market is growing steadily at a rate of over 1\% per year (missing citation); (missing citation). For instance, the global consumption of paper and board amounted to an estimated 399 million metric tons in 2020. It is expected that demand will increase steadily over the next decade, reaching approximately 461 million metric tons in 2030 (missing citation).
The vast majority of increasing demand mainly comes from the following categories of products (missing citation): personal hygiene paper products, food packaging products, corrugated packaging materials, and paper-based medicinal materials (shown in Figure \ref{paper}), which closely related to the growing e-commerce business and awareness of safety and hygiene, especially under the pandemic situation of COVID-19. 
Advanced process control and state-of-the-art process optimization techniques would provide enormous opportunities for maximized efficiency and optimized energy usage to satisfy the steadily growing need.

The pulp and paper mills aim to convert wood chips into pulp, the raw material for different types of products. In general, a pulping process can be classified into mechanical pulping, chemical pulping, and thermomechanical pulping according to the fundamental mechanism of separating wood fibers (missing citation). More than two-thirds of globally produced pulp comes from Kraft pulping which belongs to the chemical pulping process (missing citation). In the Kraft pulping process, the conversion of wood chips into pulp mainly takes place in a pressured vertical cylindrical reactor known as the pulp digester, which operates in a batch manner or as a continuous process. Due to lower space requirements and lower energy costs, continuous pulp digesters are predominantly used in industrial applications (missing citation). As illustrated in Figure \ref{digester}, the typical continuous digesters consist of three basic zones: an impregnation zone, a cook zone, and a wash zone. White liquor and pre-steamed wood chips are introduced at the top of the digester into the impregnation zone where the liquor penetrates the wet chips. After that, the two streams flow downwards into the cook zone where the most delignification reactions occur. Then, the spent liquor is withdrawn from the digester at extraction screens. At the same time, the cooked pulp moves downwards to the wash zone where the chips are washed by the counter-current flow of cold wash liquor. The cooked pulp is then removed from the bottom of the digester. In particular, the cooking degree is evaluated by the Kappa number, which denotes the residual lignin content of the pulp.

Due to the complex nature of the delignification process, numerous mathematical models for the pulp digesters have been proposed over the last few decades. Three widely used dynamic models are known as the Purdue, Gustafson, and Andersson models, which have similar conceptual bases and assumptions (missing citation); (missing citation); (missing citation). 
These kinetic models show the effect of a change of temperature on the reaction rate constants for different wood components according to the Arrhenius expression (missing citation). The main differences are the numbers of wood components, the assumption of consecutive or parallel reactions, and the assumption about how the delignification reaction takes place along the digester length. Among them, the most commonly used model is the extended Purdue model, which is followed and further developed by other researchers (missing citation); (missing citation); (missing citation).\\
When it comes to the controller design for the pulp digesters, a wide variety of control methods, such as reduced dimension control, robust control, and generic algorithms have been investigated  (missing citation); (missing citation). Since one of the main objectives of the pulp digester is to produce pulp with specific properties with minimum chemical and energy inputs, the digester is usually operated in a constrained setting(missing citation). 
 In this case, model predictive control (MPC) as a widely deployed methodology in different types of  processes is capable of handling such requirements.  Michaelsen et al. (missing citation) developed a model predictive controller for a Kamyr digester using a real-time mechanistic model compensated by an optimal state estimator and the performance demonstrated to be superior when compared with proportional-integral (PI) control in offline simulations. 
Lee and Datta  (missing citation) designed a model predictive control  system coupled with an extended Kalman filter for a batch-type pulp digester.
Along the same line of work, Wisnewski and Doyle (missing citation) analyzed the performance of linear MPC and nonlinear MPC for set-point tracking and unmeasured disturbance rejection.
Alexandridis and Sarimveis (missing citation) employed an adaptive MPC based on a radial basis function ANN model for Kappa number control of continuous pulp digesters.
H.-K. Choi and J. S.-I. Kwon (missing citation); (missing citation); (missing citation); (missing citation) developed a class of MPC for continuous pulp digester and batch pulp digester based on the proposed multi-scale model.
However, even though the aforementioned works have made a valuable contribution toward the modelling and controller design for the pulp digester, there still some aspects which did not receive much attention and/or accurate consideration. First of all, most of the works depend on the approximation of the PDE into a large-scale ordinary differential equation (ODE), which is generally prone to approximation error and is difficult to capture the spatial kinetics properties within the digester. In addition, these approaches need to perform the spatial discretization in an early lumping manner, which might induce numerical instability and/or alter the fundamental control theoretical properties, such as controllability, observability, and stabilizability (missing citation); (missing citation). Secondly, the target of the kappa number and cooking temperature setting are considerably different between softwood and hardwood pulp, the shift operations, therefore, need to be taken into account when the controller is developed. Thirdly, the co-current flow and counter-current flow in the digester need to be considered in accurate manner, mainly due to the fact that the main control manipulation is performed in the top part of the digester, while the controlled variable is physically measured at the outlet of the digester (missing citation).

Motivated by the aforementioned issues, this work considers a bilinear transformation of a continuous infinite-dimensional system to a discrete one by the application of the Cayley-Tustin time discretization approach, in which the physical characteristics (energy) and theoretical properties of considered systems are preserved (missing citation). In addition, the finite-dimensional MPC setting is extended to the distributed parameter pulp digester system to realize target tracking and achieve optimal grade transitions. The cascade zones (cook zone and wash zone) of the digester are considered as an extended distributed parameter system, which is described by 10 hyperbolic PDEs.

In this manuscript, the claimed novelty is the extension of linear MPC designs for the finite-dimensional system to the infinite-dimensional one to realize target tracking. Particularly, the pulp digester system described by a set of coupled nonlinear hyperbolic PDEs is considered, which can capture the spatiotemporal evolution of wood chips and wash liquor. In addition, the linearized continuous-time model is transformed to a linear discrete-time infinite-dimensional model, and the fundamental continuous-time properties (including stability, controllability and observability) are preserved under the Cayley--Tustin transformation. Moreover, a discrete-time observer for the system of hyperbolic PDEs is proposed, accounting for the available output measurement and the system states reconstruction. Finally, the controller design provides optimal asymptotic stabilization and target tracking realization of the system with the inclusion of output and input constraints.

The article is organized as follows: we initially present a dynamic model to describe the delignification process in the continuous pulp digester. Then, the developed model is linearized and discretized in time by utilizing the Cayley-Tustin approach. Lastly, we propose a Luenberger observer-based model predictive controller for the discrete model to realize the target tracking. The simulation results demonstrate the performance of the proposed controller design.

\section*{Problem Formulation}
In this section, the mathematical model of the pulp digester process is introduced. In particular, the equilibrium profiles are calculated to proceed with the linearization of the original nonlinear system. For the sake of simplicity, an infinite-dimensional representation is given to illustrate the linearized system. Finally, a discrete-time model is obtained utilizing the Cayley-Tustin discretization framework.

\subsection*{Model description}
The main section of a continuous pulp digester can be divided into the co-current zone and counter-current zone, which are also referred to as cook zone and wash zone respectively, as shown in Figure \ref{scheme}. Each zone can be modelled by a set of nonlinear coupled PDEs, which are derived from the first principle. The models of cook zone and wash zone are similar, but different flow direction of the liquor. Generally, each volume in the digester contains three phases, solid phase, entrapped liquor phase, and free liquor phase. In particular, the combination of solid phase and entrapped liquor phase is further referred to as the wood chip phase, and the details of the conceptual model of the mass in the digester are shown in Figure \ref{phase}. The wood chips are often assumed to consist of five components, and the entrapped liquor and free liquor contain four same components. For more details of the model, one can see the references(missing citation); (missing citation).

Based on the model proposed by Michelsen, the following assumptions are considered to obtain a simple model which describes some dominant dynamics of the process. Wood chips are assumed to be of equal size and to have a constant volume during cooking. The volume flows of wood chips and free liquor are assumed to be equal at all space positions, and the chip compaction profile is assumed to be stationary, expressed as a piecewise linear function of position. The cross-sectional area is assumed to be constant in the digester. For simplicity, the degradation of carbohydrates is assumed to follow a lignin reduction linearly and the concentrations of dissolved solids in the entrapped liquor and free liquor can be omitted  (missing citation); (missing citation); (missing citation). According to Gustafson et al. (missing citation), the carbohydrate reaction rate in the bulk delignification period is proportional to the lignin reaction rate.

Hence, from control and monitoring point of view, the following variables are important, which are considered as the states of the model:

$\rho_{s} (z, t)$: concentration of lignin of solid phase (\%)

$\rho_{e} (z, t)$: concentration of alkali of entrapped liquor phase ($g/l$)

$\rho_{f} (z, t)$: concentration of alkali of free liquor phase ($g/l$)

$T_c (z, t)$: temperature of wood chip phase (K)

$T_f (z, t)$: temperature of free liquor phase (K)

Based on the abovementioned assumption, the mathematical models of cook zone and wash zone can be simplified to the following equations:

\begin{equation}

\begin{aligned}
&0\leq z\leq L_{co}\left\{

\begin{array}
{l}
\frac{\partial \rho_{s}^{\mathcal{C}} (z, t)}{\partial t}=-v_c\frac{\partial \rho_{s}^{\mathcal{C}}  (z, t)}{\partial z}-R_{s}^{\mathcal{C}}\\ 
\frac{\partial \rho_{e}^{\mathcal{C}} (z, t)}{\partial t}=-v_c\frac{\partial  \rho_{e} ^{\mathcal{C}} (z, t)}{\partial z}+\frac{D_{EA}^{\mathcal{C}}}{\varepsilon _c}(\rho_{f}^{\mathcal{C}} (z, t)-\rho_{e}^{\mathcal{C}} (z, t))-\rho_{ODW}(b_1+0.47b_2)R_{s}^{\mathcal{C}}\\ 
\frac{\partial \rho_{f}^{\mathcal{C}}  (z, t)}{\partial t}=-v_f\frac{\partial \rho_{f}^{\mathcal{C}}  (z, t)}{\partial z}-\frac{D_{EA}^{\mathcal{C}}}{\varepsilon _f}(\rho_{f}^{\mathcal{C}} (z, t)-\rho_{e}^{\mathcal{C}} (z, t))\\ 
\frac{\partial T_c^{\mathcal{C}} (z, t)}{\partial t}=-v_c\frac{\partial T_c^{\mathcal{C}} (z, t)}{\partial z}+\left [ 1.47\rho_c \Delta H_R R_{s}^{\mathcal{C}}+U(T_f^{\mathcal{C}} (z, t)-T_c^{\mathcal{C}} (z, t))\right ] / C_{pe} \\
\frac{\partial T_f^{\mathcal{C}} (z, t)}{\partial t}=-v_f\frac{\partial T_f^{\mathcal{C}} (z, t)}{\partial z}-U(T_f^{\mathcal{C}} (z, t)-T_c^{\mathcal{C}} (z, t)) /C_{pf}+b(z)u(t)
\end{array}
\right.\\
&0\leq z \leq L_{ws}\left\{

\begin{array}
{l}
\frac{\partial \rho_{s} ^{\mathcal{W}}(z, t)}{\partial t}=-v_c\frac{\partial \rho_{s} ^{\mathcal{W}}(z, t)}{\partial z}-R_{s}^{\mathcal{W}}\\ 
\frac{\partial \rho_{e}^{\mathcal{W}}(z, t)}{\partial t}=-v_c\frac{\partial  \rho_{e}^{\mathcal{W}} (z, t)}{\partial z}+\frac{D_{EA}^{\mathcal{W}}}{\varepsilon _c}(\rho_{f}^{\mathcal{W}}(z, t)-\rho_{e}^{\mathcal{W}}(z, t))-\rho_{ODW}(b_1+0.47b_2)R_{s}^{\mathcal{W}}\\ 
\frac{\partial \rho_{f} ^{\mathcal{W}}(z, t)}{\partial t}=v_f\frac{\partial \rho_{f}^{\mathcal{W}}(z, t)}{\partial z}-\frac{D_{EA}^{\mathcal{W}}}{\varepsilon _f}(\rho_{f}^{\mathcal{W}}(z, t)-\rho_{e}^{\mathcal{W}}(z, t))\\ 
\frac{\partial T_c^{\mathcal{W}}(z, t)}{\partial t}=-v_c\frac{\partial T_c^{\mathcal{W}}(z, t)}{\partial z}+\left [  1.47\rho_c \Delta H_R R_{s}^{\mathcal{W}}+U(T_f^{\mathcal{W}}(z, t)-T_c^{\mathcal{W}}(z, t))\right ] / C_{pe} \\
\frac{\partial T_f^{\mathcal{W}}(z, t)}{\partial t}=v_f\frac{\partial T_f^{\mathcal{W}}(z, t)}{\partial z}-U(T_f^{\mathcal{W}}(z, t)-T_c^{\mathcal{W}}(z, t)) /C_{pf}
\end{array}
\right.
\label{nonpde} 

\end{aligned}

\end{equation}
where the superscript $\mathcal{C}$ and $\mathcal{W}$ indicate the cook zone and wash zone, respectively. $R_{s}^{\mathcal{C}}$ and $R_{s}^{\mathcal{W}}$ are the rate of consumption of mass of solid per chip volume in cook zone and wash zone, describing by the Arrhenius coefficients, $R_{s}^{\mathcal{C}}=A_{1}e^{-\frac{E_{1}}{T_c^{\mathcal{C}}}}\rho_{e}^{\mathcal{C}}(\rho_{s}^{\mathcal{C}}-\rho_{s}^{0 })$, $R_{s}^{\mathcal{W}}=A_{1}e^{-\frac{E_{1}}{T_c^{\mathcal{W}}}}\rho_{e}^{\mathcal{W}}(\rho_{s}^{\mathcal{W}}-\rho_{s}^{0 })$. $\rho_{s}^{0 }$ denotes the non-reactive lignin in wood. The mass diffusion rate $D_{EA}$ usually take the form, $D_{EA}^{\mathcal{C}}=0.16\sqrt{T_c^{\mathcal{C}}}e^{\frac{-2452.4}{T_c^{\mathcal{C}}}}(-2.0\rho_{s}^{\mathcal{C}}+0.13(\frac{\rho_{e}^{\mathcal{C}}}{32.0})^{0.55}+0.58)$ and $D_{EA}^{\mathcal{W}}=0.16\sqrt{T_c^{\mathcal{W}}}e^{\frac{-2452.4}{T_c^{\mathcal{W}}}}(-2.0\rho_{s}^{\mathcal{W}}+0.13(\frac{\rho_{e}^{\mathcal{W}}}{32.0})^{0.55}+0.58)$. $\varepsilon_c$ is the chip compaction, which increases from the entry through the cook zone, reaching a maximum at the main extraction, $\varepsilon_c(z)=\varepsilon _{10}+\varepsilon _{11}z$, and $\varepsilon_f(z)=1-\varepsilon_c(z)$. $U$ denotes the heat-transfer coefficient, and the heat capacities of the entrapped and free liquor phases ($C_{pe}$ and $C_{pf}$) are determined by the mixing rules based on weighted averages (missing citation). 
In this case, the input is the temperature of free liquor and the distribution is described by the actuation distribution function $b(z)=1_{[z_a-\varepsilon, z_a+\varepsilon]}$. 
The controlled output is defined as the concentration of lignin of solid phase at the outlet in wash zone, and the measurements of each state at the outlet are assumed available which are taken as the measured outputs.

The boundary conditions for the cook zone are given at $z=0$ (see Figure \ref{scheme}):

\begin{align}
\rho_{s}^{\mathcal{C}} (0,t)&=\rho_{s_0}, \rho_{e}^{\mathcal{C}} (0,t)=\rho_{e_0}, \rho_{f}^{\mathcal{C}} (0,t)=\rho_{f_0}, 
T_c^{\mathcal{C}} (0, t)=T_{c0}, T_f^{\mathcal{C}} (0, t)=T_{f0}

\end{align}
and for the wash zone at  $z=0$ and $z=L_{ws}$:

\begin{small}

\begin{align}
\rho_{s}^{\mathcal{W}}(0,t)&=\rho_{s}^{\mathcal{C}}(L_{co},t),  \rho_{e}^{\mathcal{W}}(0,t)=\rho_{e}^{\mathcal{C}}(L_{co},t),
\rho_{f}^{\mathcal{W}}(L_{ws}, t)=\rho_{f_0}, 
T_c^{\mathcal{W}}(0, t)&=T_c^{\mathcal{C}}(L_{co}, t),T_f^{wash}(L_{ws}, t)=T_{f0}
\label{bcwash} 

\end{align}

\end{small}

\subsection*{Model Linearization}
In order to reduce the mathematical complexity of the nonlinear process model, we make the linearization around the system equilibrium point or spatially nonuniform steady-state by setting the time derivative terms to be zero. The steady-states are solved numerically using the finite difference method. As a result, the corresponding steady-states profiles are illustrated in Figure \ref{ss}.
 Additionally, we assume that $L_{co}$ and $L_{ws}$ are the same, which denote as $L$ for notation simplicity, and introduce the following notations to define the local dynamics from the states of the cook zone to the states of the wash zone:

\begin{equation}

\begin{aligned}
x^{(\mathcal{C})}(z,t)&=[\rho_{s}^{\mathcal{C}}(z,t); \rho_{e}^{\mathcal{C}}(z,t); \rho_{f}^{\mathcal{C}}(z,t); T_c^{\mathcal{C}}(z,t); T_f^{\mathcal{C}}(z,t)]\\
x^{(\mathcal{W})}(z,t)&=[\rho_{s}^{\mathcal{W}}(z,t); \rho_{e}^{\mathcal{W}}(z,t); \rho_{f}^{\mathcal{W}}(z,t); T_c^{\mathcal{W}}(z,t); T_f^{\mathcal{W}}(z,t)]
 
\end{aligned}

\end{equation}

The above nonlinear model Equation (\ref{nonpde}) can be linearized applying the Taylor series expansion. The linearized model is then obtained as:

\begin{equation}

\begin{aligned}
\dot{x}^{(\mathcal{C})}(z,t)&=\mathcal{A}_1x^{(\mathcal{C})}(z,t)+\mathcal{B}_1u(t)\\
\dot{x}^{(\mathcal{W})}(z,t)&=\mathcal{A}_2x^{(\mathcal{W})}(z,t)\\
y_c(t)&=\mathcal{C}_2x^{(\mathcal{W})}(z,t)\\
y_{m}(t)&=

\begin{bmatrix}
\mathcal{C}_{m1}x^{(\mathcal{C})}(z,t)\\
\mathcal{C}_{m2}x^{(\mathcal{W})}(z,t)

\end{bmatrix}
\label{sspace} 

\end{aligned}

\end{equation}
where $z_1\in [0,L]$ and $z_2\in [0,L]$. The operator $\mathcal{A}_1(\cdot)=(V_1\frac{\partial }{\partial z}+\psi_1(z))(\cdot)$ and $\mathcal{A}_2(\cdot)=(V_2\frac{\partial }{\partial z}+\psi_2(z))(\cdot)$, where $V_1=\text{diag}\{-v_c, -v_c, -v_f, -v_c, -v_f\}$ and $V_2=\text{diag}\{-v_c, -v_c, v_f, -v_c, v_f\}$ and

\begin{small}

\begin{align}
\psi_1(z)&=

\begin{bmatrix}
J_{11}(z) &  J_{12}(z) & 0 & J_{14}(z) &0 \\ 
J_{21}(z) & J_{22}(z) & J_{23}(z) & J_{24}(z) & 0  \\ 
 J_{31}(z)&  J_{32}(z) & J_{33}(z) & J_{34}(z) &  0 \\ 
  J_{41}(z)&  J_{42}(z) & 0 & J_{44}(z) & J_{45}(z) \\ 
0& 0 & J_{53}(z) & J_{54}(z) & J_{55}(z) 
\end{bmatrix}
, \psi_2(z)=

\begin{bmatrix}
\bar{J}_{11}(z) &  \bar{J}_{12}(z) & 0 & \bar{J}_{14}(z) &0 \\ 
\bar{J}_{21}(z) & \bar{J}_{22}(z) & \bar{J}_{23}(z) & \bar{J}_{24}(z) & 0  \\ 
 \bar{J}_{31}(z)&  \bar{J}_{32}(z) & \bar{J}_{33}(z) & \bar{J}_{34}(z) &  0 \\ 
  \bar{J}_{41}(z)&  \bar{J}_{42}(z) & 0 & \bar{J}_{44}(z) & \bar{J}_{45}(z) \\ 
0& 0 & \bar{J}_{53}(z) & \bar{J}_{54}(z) & \bar{J}_{55}(z) 
\end{bmatrix}
\nonumber

\end{align}

\end{small}
where the $J_{i,j}(z)$ and $\bar{J}_{i,j}(z)$($i,j=1, 2,..., 5$) are the nonlinear term evaluated at steady-state by ignoring the 2nd order and higher order terms, such as $J_{1,1}=-A_1\rho_{e,ss}e^{-\frac{E_1}{T_{c,ss}}}$, $J_{1,2}=A_1(\rho_{s1}^{0}-\rho_{s1,ss})e^{-\frac{E_1}{T_{c,ss}}}$.
The input operator $\mathcal{B}_1$ is defined as a bounded operator $\mathcal{B}_1$ =[ 0; 0; 0; 0; $b(z)$] and $b(z)=\frac{1}{2\bar{z} _{L}}1_{[z _{L}-\bar{z} _{L},z _{L}+\bar{z} _{L}]}(z)$. The controlled output operator $\mathcal{C}_2$ is determined as $\mathcal{C}_2(\cdot) =\text{diag}\{\int_{0}^{L}\delta (z -L) (\cdot )d\eta, \!0, \!0, \!0, \!0\}$, \! \!and the measured output operators are defined as $\mathcal{C}_{m1} \! \!= \! \!\text{diag}\{\int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta\}$ and $\mathcal{C}_{m2}\!=\!\text{diag}\{\int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z ) (\cdot )d\eta, \int_{0}^{L}\delta (z -L) (\cdot )d\eta, \int_{0}^{L}\delta (z ) (\cdot )d\eta\}$, respectively.

In this form, the domains of operator $\mathcal{A}_1$ and $\mathcal{A}_2$ are $\mathcal{D}(\mathcal{A}_1) \!= \! \{\left [ \phi _1 (z); \phi _2 (z); \phi _3 (z); \phi _4 (z); \phi _5 (z)\right ] \!\in L_{2}(0,1)^5|\phi _i (z)\text{ are abs. cont.,} \text{ and }  \phi _{i}(0)\!=\!0, i=1,2,...,5 \} $, and $\mathcal{D}(\mathcal{A}_2)\!=\! \{ [ \phi _6 (z); \phi _7 (z); \phi _8 (z); \phi _9 (z);\\ \phi _{10} (z) ]\in L_{2}(0,1)^5|\phi _j (z)\text{ are abs. cont.,}  j=6,7,...,10 \text{ and } \phi _{j}(L)=0, j=8,10\} $, respectively.

Mathematically, the dynamics of the above cascade system (cook zone and wash zone) can be described by the following extended system:

\begin{equation}

\begin{aligned}
\dot{x}(z,t)&=\mathcal{A}x(z,t)+\mathcal{B}u(t)\\
y_c(t)&=\mathcal{C}_cx(z,t) \\
y_m(t)&=\mathcal{C}_mx(z,t)
\label{exsys} 

\end{aligned}

\end{equation}
with boundary conditions given below:

\begin{equation}

\begin{aligned}
x_{1} (0, t)&=0,
x_{2} (0, t)=0,
x_{3} (0, t)=0, 
x_{4} (0, t)=0,
x_{5} (0, t)=0\\
x_{6} (0, t)&=x_{1} (L, t),
x_{7} (0, t)=x_{2} (L, t),
x_{8} (L, t)=0, 
x_{9} (0, t)=x_{4} (L, t),
x_{10} (L, t)=0
\label{newbc} 

\end{aligned}

\end{equation}
where the extended state $x(z,t)=[x^{(\mathcal{C})}(z,t);x^{(\mathcal{W})}(z,t)]$ and the corresponding extended operators are given as follows:

\begin{equation}

\begin{aligned}
\mathcal{A}(\cdot)&=V\frac{\partial }{\partial z}(\cdot)+\psi(z)(\cdot), V=\text{bdiag}\{V_1,V_2\}, \psi (z)=\text{bdiag}\{\psi_1,\psi_2\}\\
\mathcal{B}&=[\mathcal{B}_1;\textbf{0} ], \mathcal{C}_c=[\textbf{0} ,\mathcal{C}_2], 
\mathcal{C}_m=\text{bdiag}\{\mathcal{C}_{m1},\mathcal{C}_{m2}\}

\end{aligned}

\end{equation}

The state $x(z,t) \in \mathcal{X}$, with $\mathcal{X}={{L}}^2((0,L),{\mathbb{R}}^{10})$ being defined as a real separable Hilbert space with inner product $\langle\cdot, \cdot\rangle$. The input $u(t) \in {{L}}_{loc}^2([0,\infty),\mathcal{U})$ and output $y(t) \in {{L}}_{loc}^2([0,\infty),\mathcal{Y}) $, where $\mathcal{U}$ and $\mathcal{Y}$ are real separable Hilbert spaces.
Based on the coupling conditions Equation (\ref{bcwash}), we have the domain of $\mathcal{A}$, which is $\mathcal{D}(\mathcal{A})=\mathcal{D}(\mathcal{A}_1)\bigoplus \mathcal{D}(\mathcal{A}_2) = \{[\phi _1;\phi _2;...;\phi _{10}]\in L_{2}(0,L)^{10}|[\phi _1;\phi _2;...;\phi _{5}] \in \mathcal{D}(\mathcal{A}_1), [\phi _6;\phi _7;...;\phi _{10}] \in \mathcal{D}(\mathcal{A}_2), \text{ and } \phi _{1}(L)=\phi _{6}(0),\phi _{2}(L)=\phi _{7}(0),\phi _{4}(L)=\phi _{9}(0) \} $ (missing citation). 
The adjoint operator $\mathcal{A}^{*}$ is easily found using the inner product formula, $\langle\mathcal{A}\varphi_i, \phi_i \rangle = \langle\varphi_i, \mathcal{A}^*\phi_i \rangle, i=1,2,...,10$, and is:

\begin{align}
\mathcal{A}^{*}(\cdot)=-V\frac{\partial {(\cdot)}}{\partial z }+\psi^{*}(z)(\cdot)

\end{align}
with its domain defined as $ \mathcal{D}(\mathcal{A}^*)=\{[\phi _1;\phi _2;...;\phi _{10}]\in L_{2}(0,1)^{10},$
$ \phi_i (z ) \text{ is absolutely continuous, } $\\$\frac{d \phi_i}{d z}\in L_{2}(0,1), \!\text{ with } i\!=\!1,2,...,10, \!\text{ and }\phi _{j}(L)\!=\!0, j=3,5,6,7,9,\phi _{m}(0)\!=\!0, m=\!8,10, \phi _{1}(L)\!=\!\frac{v_6}{v_1}\phi _{6}(0), \phi _{2}(L)=\frac{v_7}{v_2}\phi _{7}(0), \phi _{4}(L)=\frac{v_9}{v_4}\phi _{9}(0)\}$.

\section*{Model Discretization}
In this section, the Cayley-Tustin discretization framework is considered and applied to the linearized digester system without any spatial approximation induced. By the use of Cayley--Tustin transformation, a discrete-time state-space model for describing the extended digester system is established and realized by determining the resolvent operator, which is amenable to the discrete observer and controller designs. Meanwhile,  the essential properties of the continuous-time system stay invariant under this transformation, such as conservative (energy preserving)(missing citation), stability(missing citation), observability(missing citation); (missing citation), controllability.

\subsection*{Cayley-Tustin Time Discretization Framework}
The above linear system in Equation (\ref{exsys}) is considered. For a given time discretization $h>0$, and for $j\geq 1$ the Cayley-Tustin discretization is given by:

\begin{equation}

\begin{aligned}
\frac{{x(jh) - x((j - 1)h)}}{h} &\approx \mathcal{A}\frac{{x(jh) + x((j - 1)h)}}{2} + \mathcal{B}u(jh)\\
y_c(jh) &\approx \mathcal{C}_c\frac{{x(jh) + x((j - 1)h)}}{2} \\
y_m(jh) &\approx \mathcal{C}_m\frac{{x(jh) + x((j - 1)h)}}{2} 

\end{aligned}
\label{cayley} 

\end{equation}
with $x(0) = {x_0}$, where the spatial dependence of $x$ is omitted for brevity. Then, let $\frac{u_j^{(h)}}{\sqrt{h}}$ be an approximation of $u(jh)$ by the mean value within a given sampling time, $\frac{u_j^{(h)}}{\sqrt{h}}=\frac{1}{h}\int_{(j-1)h}^{jh}u(t)dt$. It has been shown in (missing citation) that $\frac{u_j^{(h)}}{\sqrt{h}}$ converges to $u(jh)$ as $h \rightarrow 0$ in several different ways, similar for $Y(jh)$. Further, rewriting Equation (\ref{cayley}) gives the discrete time dynamics Equation (\ref{tustin}). It is frequently called Tustin discretization in the engineering literature, which is proposed in 1940s by Tustin and referred as Tustin transform in digital and sample-data control literature (missing citation).

\begin{equation}

\begin{aligned}
\frac{{x_j- x_{j-1}}}{h} &{\approx \mathcal{A}\frac{{x_j+ x_{j-1}}}{2} + \mathcal{B}\frac{u_j}{\sqrt{h}}}\\
{\frac{y_{cj}}{\sqrt{h}} }&{\approx \mathcal{C}_c\frac{{x_j + x_{j-1}}}{2} }\\
{\frac{y_{mj}}{\sqrt{h}} }&{\approx \mathcal{C}_m\frac{{x_j + x_{j-1}}}{2} }

\end{aligned}
\label{tustin} 

\end{equation}
Through some basic computations, the following infinite-dimensional discrete-time state space model is obtained:

\begin{equation}

\begin{split}
x_j &= {\mathcal{A}_d}x_{j-1}+ {\mathcal{B}_d}u_j \\
y_{cj} &= {\mathcal{C}_{cd}}x_{j-1} + {\mathcal{D}_{cd}}u_j\\
y_{mj}&= {\mathcal{C}_{md}}x_{j-1}+ {\mathcal{D}_{md}}u_j

\end{split}

\end{equation}
where $\mathcal{A}_d$, $\mathcal{B}_d$, $\mathcal{C}_d$, $\mathcal{D}_d$, $ \mathcal{C}_{cd}$ and  $\mathcal{C}_{md}$ are the discrete-time spatial operators and we denote:

\begin{align}

\begin{pmatrix}
\!\!\mathcal{A}_d &  \mathcal{B}_d\\ 
 \mathcal{C}_{cd}& \mathcal{D}_{cd} \\
  \mathcal{C}_{md}& \mathcal{D}_{md}
\end{pmatrix}
\!\!=

\begin{pmatrix}
-I+2\delta \mathcal{R}(z,\delta ) & \sqrt {2\delta}\mathcal{R}(z,\delta ) \mathcal{B}\\ 
\sqrt {2\delta}\mathcal{C}_c\mathcal{R}(z,\delta ) & \mathcal{G}_c(\delta )\\
\sqrt {2\delta}\mathcal{C}_m\mathcal{R}(z,\delta ) & \mathcal{G}_m(\delta )
\end{pmatrix}
\label{operators} 

\end{align}
where $\delta=2/h$ and the resolvent is $\mathcal{R}(z,\delta)=(\delta -\mathcal{A})^{-1}$. Clearly, one must satisfy $\delta \in \rho(\mathcal{A})$ so that the resolvent operator is well-defined. In particular, $\mathcal{G}_c(\delta )$ denotes the transfer function from input to controlled output $\mathcal{G}_c(\delta )=\mathcal{C}_c(\delta  - \mathcal{A})^{ - 1}\mathcal{B}$, and $\mathcal{G}_m(\delta )$ denotes the transfer function from input to measured outputs $\mathcal{G}_m(\delta )=\mathcal{C}_m(\delta  - \mathcal{A})^{ - 1}\mathcal{B}$.
The unbounded operators $\mathcal{A} $ of the continuous-time system are mapped into bounded operators $\mathcal{A}_d$ in the discrete-time counterpart through Cayley transform, which describe by the resolvent operator with $\mathcal{A}_d(\cdot )=[\delta  -\mathcal{A}]^{-1}[\delta  +\mathcal{A}](\cdot )=-I(\cdot )+2\delta[\delta  -\mathcal{A}]^{-1}(\cdot )=-I(\cdot )+2\delta\mathcal{R} (\delta,\mathcal{A})(\cdot )$. In addition, it has been demonstrated that the controllability and stability are invariant under this transformation (missing citation).

\subsection*{Resolvent operator}
The resolvent operator can be obtained by utilizing the Laplace transform. Under the zero input condition, we can have the following expression,

\begin{align}
\frac{\partial x(z,s)}{\partial z}=V^{-1}(sI-\psi)x(z,s)-V^{-1}x(z,0)

\end{align}
where $\psi$ is taken as the average of $\psi(z)$. By solving the above ODE, one obtains:

\begin{align}
x(z,s)=e^{\{V^{-1}(sI-\psi)z\}}x(0,s) -\int_{0}^{z}e^{\{V^{-1}(sI-\psi)(z-\eta)\}}V^{-1}x(\eta ,0)d\eta
\label{lap} 

\end{align}
Since $\psi$ is the block diagonal matrix and $V$ is already in diagonal form, the matrix $F$ is finally a block diagonal matrix, denoted as follows:

\begin{align}
F=V^{-1}(sI-\psi)=

\begin{bmatrix}
[F_1]_{5\times 5} & \mathbf{0}_{5\times 5}\\ 
\mathbf{0}_{5\times 5} & [F_2]_{5\times 5}

\end{bmatrix}

\end{align}
Then $e^{Fz}=\text{diag}(e^{F_1z}, e^{F_2z})$ and we denote $[M_{ij}(z,s)]_{10\times 10}=e^{Fz }$ for simplicity.

In order to obtain the unknown boundary conditions, one needs to evaluate the boundary conditions in Equation (\ref{newbc}) as follows:

(1) At $z = 0$, one can plug $x_{1} (0, t)=0, x_{2} (0, t)=0, x_{3} (0, t)=0,  x_{4} (0, t)=0, x_{5} (0, t)=0$ into Equation (\ref{lap}) which leads to $M_{ij}(0,s)=0, i=1,2,...5, j=6,7,...10$.

(2) At $z = L$, one can firstly substitute $x_{6} (0, s)=x_{1} (L, t), x_{7} (0, s)=x_{2} (L, t),x_{9} (0, s)=x_{4} (L, t)$ into Equation (\ref{lap}) which yields

\begin{equation}

\begin{aligned}
x_{6} (0, s)=-\sum_{k=1}^{10}\int_{0}^{L}M_{1k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \\
x_{7} (0, s)=-\sum_{k=1}^{10}\int_{0}^{L}M_{2k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \\
x_{9} (0, s)=-\sum_{k=1}^{10}\int_{0}^{L}M_{4k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta

\end{aligned}

\end{equation}
where $V_v$ denotes $V^{-1}$.
Then, substituting the boundary conditions $x_{8} (L, s)=0,x_{10} (L, s)=0$, one can have

\begin{small}

\begin{align}
&

\begin{bmatrix}
 M_{88}&  M_{810}\\ 
  M_{108}&  M_{1010}
\end{bmatrix}

\begin{bmatrix}
x_8(0,s)\\ 
x_{10}(0,s)
\end{bmatrix}
\nonumber\\ 
&=

\begin{bmatrix}
\sum_{k=1}^{10}\int_{0}^{L}M_{8k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta-M_{86}x_6(0,s)-M_{87}x_7(0,s)-M_{89}x_9(0,s)\\ 
\sum_{k=1}^{10}\int_{0}^{L}M_{10k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta-M_{106}x_6(0,s)-M_{107}x_7(0,s)-M_{109}x_9(0,s)

\end{bmatrix}

\end{align}

\end{small}
By substituting $x_{6} (0, s)$, $x_{7} (0, s)$, $x_{9} (0, s)$ and solving the above equations, we can obtain the expressions of $x_{8} (0, s)$ and $x_{10} (0, s)$ as follows,

\begin{align}
x_8(0,s)=&G_{11}(\sum_{k=1}^{10}\int_{0}^{L}M_{8k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta)+G_{12}(\sum_{k=1}^{10}\int_{0}^{L}M_{10k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta) \nonumber\\ 
&+(G_{11}M_{86}+G_{12}M_{106})\sum_{k=1}^{10}\int_{0}^{L}M_{1k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \nonumber \\
&+(G_{11}M_{87}+G_{12}M_{107})\sum_{k=1}^{10}\int_{0}^{L}M_{2k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \nonumber\\
&+(G_{11}M_{89}+G_{12}M_{109})\sum_{k=1}^{10}\int_{0}^{L}M_{4k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \\
x_{10}(0,s)=&G_{21}(\sum_{k=1}^{10}\int_{0}^{L}M_{8k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta)+G_{22}(\sum_{k=1}^{10}\int_{0}^{L}M_{10k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta) \nonumber\\ 
&+(G_{21}M_{86}+G_{22}M_{106})\sum_{k=1}^{10}\int_{0}^{L}M_{1k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \nonumber\\
&+(G_{21}M_{87}+G_{22}M_{107})\sum_{k=1}^{10}\int_{0}^{L}M_{2k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta \nonumber\\
&+(G_{21}M_{89}+G_{22}M_{109})\sum_{k=1}^{10}\int_{0}^{L}M_{4k}(L-\eta,s )V_{vkk}x_k(\eta ,0)d\eta

\end{align}
where

\begin{align}
G=

\begin{bmatrix}
 M_{88}&  M_{810}\\ 
  M_{108}&  M_{1010}
\end{bmatrix}
^{-1}

\end{align}

Therefore, the resolvent operator is determined as following form by substituting the boundary conditions (Equation (\ref{newbc})):

\begin{align}
x(z,s)=[\mathcal{R}_{ij}(z,s)]_{10\times 10}x(\eta,0)

\end{align}
where

\begin{align}
\mathcal{R}_{ij}=\left\{

\begin{matrix}
\mathcal{R}c_{ij}, i\in [1,5], j\in [1,10]\\ 
\mathcal{R}w_{ij}, i\in [6,10], j\in [1,10]
\end{matrix}
\right.

\end{align}
where

\begin{align}
\mathcal{R}c_{ij}=-\int_{0}^{z}M_{ij}(z-\eta ,s)V_{vjj}(\cdot )d\eta, i\in [1,5], j\in [1,10] 

\end{align}

\begin{small}

\begin{equation}

\begin{aligned}
\mathcal{R}w_{ij}&=(M_{i8}(z,s)G_{11}+M_{i10}(z,s)G_{21})\int_{0}^{L}M_{8j}(L-\eta,s )V_{vjj}(\cdot )d\eta \nonumber\\
&+(M_{i8}(z,s)G_{12}+M_{i10}(z,s)G_{22})\int_{0}^{L}M_{10j}(L-\eta,s )V_{vjj}(\cdot )d\eta \nonumber\\
&+(M_{i8}(z,s)G_{11}M_{86}(L,s)+M_{i8}(z,s)G_{12}M_{106}(L,s)+M_{i10}(z,s)G_{21}M_{86}(L,s) \nonumber\\
&+M_{i10}(z,s)G_{22}M_{106}(L,s)-M_{i6}(z,s))\int_{0}^{L}M_{1j}(L-\eta,s )V_{vjj}(\cdot )d\eta
  \nonumber\\
&+(M_{i8}(z,s)G_{11}M_{87}(L,s)\!+\!M_{i8}(z,s)G_{12}M_{107}(L,s)\!+\!M_{i10}(z,s)G_{21}M_{87}(L,s)\!+\!M_{i10}(z,s)G_{22}M_{107}(L,s) \nonumber\\
&-M_{i7}(z,s))\int_{0}^{L}M_{2j}(L-\eta,s )V_{vjj}(\cdot )d\eta+(M_{i8}(z,s)G_{11}M_{89}(L,s)+M_{i8}(z,s)G_{12}M_{109}(L,s) \nonumber\\
&+M_{i10}(z,s)G_{21}M_{89}(L,s)+M_{i10}(z,s)G_{22}M_{109}(L,s)-M_{i9}(z,s))\int_{0}^{L}M_{4j}(L-\eta,s )V_{vjj}(\cdot )d\eta \nonumber\\
&-\int_{0}^{z}M_{ij}(z-\eta ,s)V_{vjj}(\cdot )d\eta, i\in [6,10], j\in [1,10]

\end{aligned}

\end{equation}

\end{small}
The discrete-time operators in Equation (\ref{operators}) can be solved by straightforwardly substituting the above resolvent operators. Afterwards, the discrete-time linear model is obtained:

\begin{equation}

\begin{split}
x(z ,k) &= {\mathcal{A}_d}x(z ,k - 1) + {\mathcal{B}_d}u(k)\\
y_c(k) &= {\mathcal{C}_{cd}}x(z ,k - 1) + {\mathcal{D}_{cd}}u(k)\\
y_m(k) &= {\mathcal{C}_{md}}x(z ,k - 1) + {\mathcal{D}_{md}}u(k)

\end{split}
\label{discretetimemodel} 

\end{equation}
with the boundary conditions (Equation (\ref{newbc})).

\section*{Observer-Based MPC Design}
An observer-based model predictive controller is designed for the discrete-time pulp digester system. In particular, a discrete Luenberger observer is designed first to reconstruct the states based on the available real-time measurements. The Luenberger observer is one of the practical and easy-to-realize observer, which is further considered in a discrete setting controller realization.
 The constrained optimal controller design for the finite-dimensional system theory is extended and deployed for the infinite-dimensional digester system. The overall closed-loop operation of the digester process is schematically presented in Figure \ref{framework}.

\subsection*{Discrete Luenberger observer design}
Firstly, let us recall the linearized discrete-time model:

\begin{equation}

\begin{split}
x(z ,k) &= {\mathcal{A}_d}x(z ,k - 1) + {\mathcal{B}_d}u(k)\\
y_m(k) &= {\mathcal{C}_{md}}x(z ,k - 1) + {\mathcal{D}_{md}}u(k)

\end{split}
\label{dtmodel} 

\end{equation}
The discrete Luenberger observer is presented by the following standard form:

\begin{equation}

\begin{split}
\hat{x}(\zeta ,k) &= {\mathcal{A}_d}\hat{x}(\zeta ,k - 1) + {\mathcal{B}_d}u(k) +L_d(y_m(k)-\hat{y}_m(k)) \\
\hat{y}_m(k)&={\mathcal{C}_{md}}\hat{x}(\zeta ,k - 1) + {\mathcal{D}_{md}}u(k)

\end{split}
\label{disluen} 

\end{equation}
where the reconstructed state $\hat{x}(\zeta ,k)$ is defined as a copy of
the system dynamics and $L_d$ is the discrete observer gain to be designed. Stability of the observer implies that the state estimation error, $e_k=x(\zeta, k)-\hat{x}(\zeta, k)$, converges to zero within a certain time. The error dynamic equation is shown as follows:

\begin{align}
e_k=(\mathcal{A}_d-L_d\mathcal{C}_d)e_{k-1}=\tilde{\mathcal{A}}_de_{k-1}
\label{errordy} 

\end{align}
To guarantee the operator $\tilde{\mathcal{A}}_d$ in the state
estimation error dynamics given by Equation (\ref{errordy}) is stable, the design objective is to determine the appropriate spatially varying gain $L_d$.
By Curtain and Zwart (missing citation), it can be shown that the operator $\tilde{\mathcal{A}}_d$ is power stable if and only if there exists a non-negative self-adjoint operator $Q_d$ such that

\begin{align}
\tilde{\mathcal{A}}_dQ_d\tilde{\mathcal{A}}_d^*-Q_d=- M_d, \text{on }    \mathcal{X}

\end{align}
where $M_d$ is a positive definite design parameter.

Then, let us assume that the pair $(\mathcal{A}, \mathcal{C})$ is exponential detectable, then if there exists a nonnegative self-adjoint operator $Q_d$ which is the solution of the following operator Riccati equation (missing citation); (missing citation):

\begin{align}
\mathcal{A}_dQ_d\mathcal{A}_d^*-Q_d-L_d(P+2I)L_d^*+M_d=0, \text{on }    \mathcal{X}
\label{disricca} 

\end{align}
where

\begin{align}
L_d=T(P+I)^{-1},
T=\mathcal{A}_dQ_d\mathcal{C}_d^*,
P=\mathcal{C}_dQ_d\mathcal{C}_d^*  \nonumber

\end{align}
The observer gain $L_d= \mathcal{A}_dQ_d\mathcal{C}_d^*(\mathcal{C}_dQ\mathcal{C}_d^*+I)^{-1}$ is a strongly stabilizing gain which guarantees the power stability of $\tilde{\mathcal{A}}_d=\mathcal{A}_d-L_d\mathcal{C}_d$. To solve the algebraic Riccati Equation (\ref{disricca}), one can utilize the numerical iteration methods, such as Newton-Kleinman iteration method (missing citation), and the detailed procedures were provided (missing citation).

\subsection*{MPC design for target tracking}
The ultimate objective of a pulp and paper mill is to ensure the specified quality of the end products while meeting the production targets and minimizing the operational costs. Specifically, when it comes to practical applications, the most characteristic point is that the species of feed wood chips are switched frequently, softwood and hardwood chips, to supply the required amount of pulp according to the production schedule of the paper machines (missing citation). Accordingly, the desired product quality needs to change in the middle of the operation, i.e., the shift operations. To realize it, the MPC is developed for the infinite-dimensional setting, emerging from the finite-dimensional linear time-invariant systems,
see Rawlings et al. (missing citation).

In this case, we consider that the system output is required to track a nonzero target vector, $y_t$, then state and input vectors, $x_t$ and $u_t$, are required which bring the system to $y_t$ at steady-state (missing citation); (missing citation); (missing citation). The state and input target can be computed by solving the following quadratic program.

\begin{equation}

\begin{split}
\mathop {\min }\limits_{{x_t, u_t}} \;  & { (u_t - \bar{u}) ^TR_t(u_t - \bar{u}) }\\
{\rm{s.t. \;\;}}&

\begin{bmatrix}
I-A_d & -B_d\\ 
C_d & D_d
\end{bmatrix}

\begin{bmatrix}
x_t\\ 
u_t
\end{bmatrix}
=

\begin{bmatrix}
0\\ 
y_t
\end{bmatrix}
 \\
&{u_{\min }} \le u_t \le {u_{\max }} \\
&{y_{\min }} \le y_t \le {y_{\max }} 

\end{split}
\label{stateinput} 

\end{equation}

In this quadratic program, $\bar{u}$ is the set point for the manipulated variables and $R_t$, is symmetric positive definite. Notice that, often, the input set point is not specified and it can be assumed zero in order to use (missing citation). The equality constraints guarantee a steady-state solution and offset free tracking of the target vector (missing citation). $u_{min}$, $u_{max}$, $y_{min}$, $y_{max}$ are the input and output constraints, respectively.

Then, the the following quadratic objective function is used for the regulator to track a nonzero target vector.

\begin{equation}

\begin{split}
\mathop {\min }\limits_{{u^N}} \;  &\sum\limits_{j = 0}^{N - 1}\;  { {\langle ( y_{k + j|k}-y_t ) ,Q(y_{k + j|k} -y_t) \rangle}  + \langle (u_{k + j|k}-u_t), R(u_{k + j|k}-u_t)\rangle} \\
& + \langle (x_{k + N-1|k}-x_t), \bar Q(x_{k + N-1|k}-x_t)   \rangle  \\
{\rm{s.t. \;\;}}&\text{Lunberger observer, Equation }(\ref{disluen})\\
&{u_{\min }} \le u_{k + j|k} \le {u_{\max }} \\
&{y_{\min }} \le y_{k + j|k} \le {y_{\max }} 

\end{split}
\label{trackmpc} 

\end{equation}
where $N$ is the prediction horizon, and $Q$, $R$ are symmetric positive semidefinite and symmetric positive definite spatial operator, respectively. $ \bar{Q}$ is the spatial operator to penalize the terminal state which depends on the stability of the given model. 
The target vector $x_t$, and $u_t$ are computed from the quadratic program in Equation (\ref{stateinput}). $y_{k + j|k}$ and $u_{k + j|k}$ represent the output and input variables at future time $k+j$ predicted at current time $k$, and the term $\Delta u_{k + j|k}$ denotes the change of an input vector at time $k+j$ as $\Delta u_{k + j|k} = u_{k + j|k}-u_{k + j-1|k}$. The vector $u^N$ includes the control sequence $\{u_{k|k}, u_{k + 1|k}, u_{k + 2|k}, ..., u_{k + N-1|k}\}$ and the first element $u_{k|k}$ will be injected to the plant as the future control action.

As discussed by H. Kwakernaak and R. Sivan (missing citation), using the targets computed from Equation (\ref{stateinput}), we define a shifted input $\tilde{u}_k=u_k - u_t$, a shifted state $\tilde{x}_k=\hat{x}_k - x_t$, and a shifted output $\tilde{y}_k=\hat{y}_k - y_t$ to reduce the problem to the standard form. The corresponding constraints can be translated to constraints on $\tilde{u}$ and $\tilde{y}$.

Thus the regulator optimization problem Equation (\ref{trackmpc}) becomes:

\begin{equation}

\begin{split}
\mathop {\min }\limits_{{u^N}} \;  &\sum\limits_{j = 0}^{N - 1}\;  { {\langle \tilde{y}_{k + j|k} ,Q\tilde{y}_{k + j|k} \rangle}  + \langle \tilde{u}_{k + j|k}, R\tilde{u}_{k + j|k}\rangle}  + \langle \tilde{x}_{k + j|k}, \bar Q\tilde{x}_{k + j|k}   \rangle  \\
{\rm{s.t. \;\;}}&\tilde{x}_{k + j|k} = {\mathcal{A}_d}\tilde{x}_{k + j-1|k} + {\mathcal{B}_d}\tilde{u}_{k + j|k}  \\
&\tilde{y}_{k + j|k} = {\mathcal{C}_d}\tilde{x}_{k + j-1|k} + {\mathcal{D}_d}\tilde{u}_{k + j|k} \\
&{u_{\min }-u_{t|k}} \le \tilde{u}_{k + j|k} \le {u_{\max }-u_{t|k}} \\
&{y_{\min }-y_{t|k}} \le \tilde{y}_{k + j|k} \le {y_{\max }-y_{t|k}} 

\end{split}
\label{trackmpc2} 

\end{equation}
According to the nature of transport reaction systems, the operator $ \bar{Q}$ can be determined from the positive self-adjoint solution of the following discrete-time Lyapunov equation:

\begin{align}
{\mathcal{A}^*_d}\bar Q \mathcal{A}_d  - \bar Q
=- {\mathcal{C}^*_d}Q\mathcal{C}_d
\label{ctlya} 

\end{align}
or equivalently the continuous-time Lyapunov equation (missing citation):

\begin{align}
{\mathcal{A}^*}\bar Q + \bar Q\mathcal{A} =  - {\mathcal{C}^*}Q\mathcal{C}
\label{dtlya} 

\end{align}
on the dual space of $\mathcal{X}_{-1}$. In addition, the operator $ \bar{Q}$ is the unique positive self-adjoint solution of the Lyapunov equations (Equation (\ref{ctlya}) and (\ref{dtlya})) (missing citation); (missing citation).

Before further manipulate the objective function (Equation (\ref{trackmpc2})), we introduce the following notations: $Y_k=\{y_{k+n}\}_{n=1}^{N}\in \mathcal{Y}^{N}$ and $U_k=\{u_{k+n}\}_{n=1}^{N}\in \mathcal{U}^{N}$. As a result, the straightforward algebraic manipulation of the objective function presented in Equation (\ref{trackmpc2}) leads to the following quadratic programming optimization problem:

\begin{equation}

\begin{split}
\mathop {\min }\limits_{U_k} \;& \langle U_k, HU_k \rangle + 2\langle U_k,F\tilde{x}_k  \rangle  +   \langle \tilde{x}_k,\bar Q\tilde{x}_k  \rangle  \\
{\rm{s.t. \;\;}}&{U_{\min }} \le U_k \le {U_{\max }} \\
&{Y_{\min }} \le GU_k + S\tilde{x}_k \le {Y_{\max }}

\end{split}
\label{qprog} 

\end{equation}
where $H \in \mathcal{L}( \mathcal{U}^N)$ is positive and self-adjoint, which is given by:

\begin{align}
h_{i,j}=\left\{

\begin{matrix}
\mathcal{D}_{d}^{*}Q\mathcal{D}_d+\mathcal{B}_{d}^{*}\bar{Q}\mathcal{B}_d+R & \text{for }  i=j \\
\mathcal{D}_{d}^{*}Q\mathcal{C}_d\mathcal{A}_{d}^{i-j-1}\mathcal{B}_d+\mathcal{B}_{d}^{*}\bar{Q}\mathcal{A}_{d}^{i-j}\mathcal{B}_d &  \text{for } i>j \\
h_{j,i}^{*} &  \text{for } i<j
\end{matrix}
\right.

\end{align}
and $F$ is given by $F=\{\mathcal{D}_d^*Q{\mathcal{C}_d}\mathcal{A}_{d}^{k-1}+\mathcal{B}_d^*\bar Q\mathcal{A}_d^k\}_{k=1}^{N-1}$. The mateix $G$ is a lower triangular given by

\begin{align}
g_{i,j}=\left\{

\begin{matrix}
\mathcal{D}_{d} & \text{for }  i=j \\
\mathcal{C}_d\mathcal{A}_{d}^{i-j-1}\mathcal{B}_d &  \text{for } i>j \\
0 &  \text{for } i<j
\end{matrix}
\right. \nonumber

\end{align}
and $S=\{{\mathcal{C}_d}\mathcal{A}_{d}^{k-1}\}_{k=1}^N$.

The inner products in the objective function given in Equation (\ref{qprog}) are vector products as $U$ is the finite-dimensional input space, and therefore we have a finite dimensional quadratic optimization problem:

\begin{equation}

\begin{split}
\mathop {\min }\limits_{U_k} \;&U_k^THU_k  + 2U_k^TF\tilde{x}_k    \\
{\rm{s.t. \;\;}}&

\begin{bmatrix}
I\\ 
-I\\ 
G\\ 
-G
\end{bmatrix}
U_k \le

\begin{bmatrix}
U_{\max }\\
-U_{\min }\\
Y_{\max }-S\tilde{x}_k\\
-Y_{\min }+S\tilde{x}_k

\end{bmatrix}

\end{split}
\label{qprognew} 

\end{equation}
Here we neglect the term $\langle \tilde{x}_k,\bar Q\tilde{x}_k  \rangle$ as $\tilde{x}_k$ is the initial condition for step $k+1$ and cannot be affected by the control input. Thus, the optimal input trajectory ($U_k$) can be obtained as the solution of the feasible quadratic optimization problem (Equation (\ref{qprognew})) converges to zero.

\section*{Numerical Simulations}
In this section, the closed-loop performance of the proposed MPC framework is demonstrated. The temperature of free liquor flowing into the cook zone is selected as the manipulated input variable and the concentration of lignin is selected as the controlled output variable. In this case, we consider the target tracking of the output by using the proposed MPC. The resulting constrained optimization problems become quadratic programming problems which are solved using the MATLAB subroutine QuadProg. Both the control and prediction horizons are chosen to be 50 sampling periods. The sampling time is set to be 10 min and the internal spatial discretization is taken as 0.05 m.

The values of all system parameters taken in the simulations are listed in Table \ref{tab1}. For the initial conditions of the dynamic system, we consider $x_{1}(z, 0)=0.067z, x_{2}(z, 0)=0.1646z, x_{3}(z, 0)=0.2377z, x_{4}(z, 0)=1.7073z, x_{5}(z, 0)=0.8661z, x_{6}(z, 0)=0.0336+2.87\times10^{-4}z, x_{7}(z, 0)=0.8232+0.0095z, x_{8}(z, 0)=0.0012z, x_{9}(z, 0)=8.5363+0.1767z$, and  $x_{10}(z, 0)=-0.086z$. In addition, $\bar{\zeta} _{L}=0.15$ is chosen for the input operator. 
 The selected parameter values for MPC implementation are listed in Table \ref{tab2} .
 
As illustrated in Figure \ref{x}, the open-loop states converge to their corresponding steady-states rapidly which indicates the original plant is intrinsically stable. The figures in the left column represent the evolution of the states in the cook zone, and the right column denotes the evolution of the states in the wash zone.
Figure \ref{x}(a) and \ref{x}(f) shows how the concentration of lignin decreases smoothly down toward the end of the cooking zone where the reactions
are stopped (or quenched) by displacement of the hot liquor with dilute wash liquor from below. Hence, no significant decrease occurs in the wash zone. 
The concentration of alkali of entrapped liquor phase in the feed flow to the cook zone is about $23$ g/l (as shown in Figure \ref{x}(b)) and then is consumed giving a decreasing profile down toward the extraction. At the bottom of the wash zone, the concentration of alkali of entrapped liquor phase is about $8$ g/l (as shown in Figure \ref{x}(g)). A similar trend occurs for the alkali of the free liquor phase, as shown in Figures \ref{x}(c) and \ref{x}(h).
Figure \ref{x}(d) and \ref{x}(e) shows the temperature profiles for the two phases (wood chip phase and free liquor phase) at cook zone, and the temperature profiles of them at wash zone shown in Figure \ref{x}(i) and \ref{x}(j). The temperature of wood chips rises due to the exothermic reactions and is also affected by the high temperature of the circulation liquor. Below the extraction screens, in the wash zone, the chips are rapidly cooled down by the wash water which has a temperature of $360$K at the inlet in the bottom.

By implementing the proposed MPC frameworks, we aim to steer the system to the desired targets without violating the physical constraints of actuators and sensors. In this case, we consider the shifted output targets, which are chosen as $y_{t1}=0.06$ and $y_{t2}=0.14$, respectively.
First of all, the steady-state target $x_t$ and input vector $u_t$ are computed from the quadratic program Equation (\ref{stateinput}). Then, based on the operator Riccati Equation (\ref{disricca}), we determine $Q_d$ in the discrete-time observer design. The estimate state $\hat{x}$ from the Luenberger observer is then utilized for the MPC design. Finally, the optimal input trajectory is computed by solving the above optimization problem outlined in Equation (\ref{qprognew}) with a receding prediction horizon strategy.
 The state evolutions of the closed-loop system under the MPC law are obtained and shown in Figure \ref{closex}. Comparing with the open-loop state profiles, the closed-loop system is able to track the target steady-states when the grade shift happens ($t=370$ min).
In addition, three pairs of states in cook zone and wash zone are successfully connected through the extend system (Equation (\ref{exsys})), such as, $\rho_{s}^{\mathcal{C}}(z,t)$ and $\rho_{s}^{\mathcal{W}}(z,t)$; $\rho_{e}^{\mathcal{C}}(z,t)$ and $\rho_{e}^{\mathcal{W}}(z,t)$; $T_c^{\mathcal{C}}(z,t)$ and $T_c^{\mathcal{W}}(z,t)$.

The open-loop output profile and the close-loop profile under the observer-based MPC law are shown in Figure \ref{tr}. Without implementing the controller, the concentration of lignin converges to its steady-state, while violating the given constraints of the system. In the closed-loop system, the concentration of lignin is able to track the targets or desired values and satisfies the requirements of the constraint simultaneously. As the targets are switched, the output can also achieve target tracking across the original steady-state.
Specifically, the target tracking above the steady-state is realized in the first period, $t\in[0,370]$, where $y_{t1}=0.06$ is considered. Similarly, the target tracking blows the steady-state is realized in the second period (i.t.,$t\in(370,740]$) when $y_{t2}=0.14$ is taken into account.
Typically, this results in higher yield, thus lowering the operating cost significantly (missing citation).
The free-liquor temperature profile, computed by the proposed model-based MPC system at each sampling time, is presented in Figure \ref{u}. 
The input trajectory corresponds to the output variables, that is, in the first 370 minutes, the optimal input variables fluctuate between 425-430 and converge to 427K when the output goes to track target 1. In the last 370 minutes, the corresponding input rises and stabilizes at 436K when target 1 switches to target 2. In these two stages of tracking, the input variables are constrained within the given bounds of actuators. The simulation studies demonstrate that the extended system is able to describe the dynamics of the original cascade system which contains a co-current zone and a counter-current zone. Furthermore, it is possible to realize the optimal control of the output in the wash zone through the operation of the input in the cook zone. The effectiveness can be demonstrated from the proposed MPC design.

\section*{Conclusions}
In this work, dynamic modelling and model predictive control design of a continuous pulp digester described by ten linearized first-order coupled hyperbolic equations was developed. The connected cook zone and wash zone of the digester were modelled as a cascade PDE system. 
By using Cayley--Tustin transformation, the linearized continuous-time infinite-dimensional model was transformed into a discrete-time infinite-dimensional model without spatial discretization and model reduction which preserves the input-output stability of the system. A Luenberger observer was designed to realize the state estimation of the system and the discrete-time Riccati equation was used to calculate the observer gain. The model predictive controller was formulated on that basis to realize target tracking and account for input and output constraints when it comes to the shift operations of the digester. The closed-loop simulation results have demonstrated that the controlled variables were able to reach to the target values and satisfy the constraints of actuators simultaneously.

\newpage

\newpage

\newpage
\listoffigures
\newpage

\newpage\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/papes/papes}
\caption{{\label{paper} The categories of paper products%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/digestes/digestes}
\caption{{\label{digester} Simplified scheme of a continuous pulp digester \protect(missing citation)%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/newschf/newschf}
\caption{{\label{scheme} The cook zone and wash zone in a digester%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/phaserepresenu/phaserepresenu}
\caption{{\label{phase} The conceptual model of the mass in a digester%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/ssopeo/ssopeo}
\caption{{\label{ss} Steady-state profiles of the digester (The solid lines denote the steady-states of the components in cook zone; the dash lines denote the steady-states of the corresponding components in wash zone.)%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/framenex/framenex}
\caption{{\label{framework} The proposed closed-loop operation framework%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/openx3/openx3}
\caption{{\label{x} The open-loop state profiles of the digester%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/closex6/closex6}
\caption{{\label{closex} The state profiles of the digester under closed-loop operation%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/trackobsvstair3/trackobsvstair3}
\caption{{\label{tr} The concentration profile of lignin under closed-loop operation%
}}
\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/tracku3/tracku3}
\caption{{\label{u} The optimal manipulated input trajectory under closed-loop operation%
}}
\end{center}
\end{figure}

\newpage
\listoftables
\newpage

\newpage\selectlanguage{english}
\begin{table}[H]
\caption {{Notations and values of parameters }} \label{tab1} 
\vspace{1.6mm}
\centering
%\begin{small}
\begin{tabular}{c c c c c c c c c }
\hline
 \hline
Process parameters& Notations & Numerical Values \\
\hline
Volumetric flow rate of chip     & $\dot{V}_c $        & 0.0267 $m^3/min      $  \\
%(1.3964)             \\ 
Volumetric flow rate of free liquor      & $\dot{V}_f  $    & 0.09 $ m^3/min         $\\              
%(2.3497)\\ 
Digester cross sectional area     &$ A    $     & 1 $ m^2      $                     \\
Concentration of non-reactive lignin &$ \rho_{s}^{0}    $     & 0.015 $ kg/m^3      $                     \\
Frequency factor of lignin reactions &$ A_1    $     & 0.09 $ m^3/kg\cdot min      $                     \\ 
Activation energy for lignin &$ E_1    $     & 38 $ kJ/mol \cdot K      $                     \\ 
Stoichiometric coefficient for lignin reactions &$  b_1   $     & 0.15                     \\ 
Stoichiometric coefficient for carbohydrates reactions &$  b_2   $     & 0.25                     \\ 
Heat capacities of the wood &$C_{ps}    $     & 1.47 $  kJ/kg\cdot K      $                     \\
Heat capacities of the liquor &$C_{pl}    $     & 4.19 $  kJ/kg\cdot K      $                     \\
Heat of reaction&$\Delta H_R    $     & -581 $  kJ/kg      $                     \\
Water density&$\rho_w    $     & 1000 $  kg/m^3      $                     \\
 \hline
  \hline
\end{tabular}
%\end{small}

\end{table}

\vspace{-1.6mm}\selectlanguage{english}
\begin{table}[H]
\caption {{Parameters for the MPC Design }} \label{tab2} 
\vspace{1.6mm}
\centering
\begin{tabular}{c c c}
\hline
 \hline
Descriptions       & Notations                   & Values            \\
\hline
Sampling time      & $h $                          & 10 $min$            \\
Prediction horizon & $N $                          & 50                \\
Input weight       & $R$                           & 0.5               \\
Output weight      & $Q$                           & 0.5               \\
Input constraints   & $[u_{min},u_{max}]$ & $[414.5K, 440.5K]$  \\
Controlled output constraints  &$ [y_{min},y_{max}{]}$ & $[0.04\%, 0.26\%]$ \\
 \hline
  \hline
\end{tabular}

\end{table}

\selectlanguage{english}
\FloatBarrier
\section*{References}\sloppy
\end{document}

