566
Ind. Eng. Chem. Res. 2004, 43, 566-576
Manufacture of Thin-Film Solar Cells: Modeling and Control of Cu(InGa)Se2 Physical Vapor Deposition onto a Moving Substrate S. Tobias Junker Department of Chemical Engineering, University of Delaware, Newark, Delaware 19716
Robert W. Birkmire Institute of Energy Conversion, University of Delaware, Newark, Delaware 19716
Francis J. Doyle, III* Department of Chemical Engineering, University of California, Santa Barbara, California 93106
Absorber layers for thin-film Cu(InGa)Se2 solar cells are produced by co-evaporation from elemental effusion sources. To meet the product specifications on thickness and composition, the vapor flow rates from the sources need to be controlled tightly. The thickness and composition of the final product are measured by X-ray fluorescence (XRF). Using the cinematic technique, a discrete dynamic model of the deposition process has been developed, and together with the XRF model, it can be cast into a Wiener model structure. The nonlinear map is globally invertible, which allows for linearization and complete decoupling of the individual control loops. Two types of controllers were tested in simulation studies: model predictive control (MPC) and internal model control (IMC). For MPC, both controllers with output bias and state estimation were developed; for IMC, both first-order plus time delay and state space models were used. Controller performance for reference tracking is compared under input, parametric, and output uncertainty. 1. Introduction Presently, most solar cell systems are based on crystalline silicon technology, which is quickly approaching its practical limit in terms of production costs. Thin-film solar cells based on Cu(InGa)Se2 absorber layers are a promising new technology allowing for economical high-volume manufacturing techniques and thus dramatically reducing costs. These technologies are in a premanufacturing development stage with largely custom-designed equipment and very complex processes. As a result, scale-up has been much more difficult than expected.1 To date, the highest efficiency solar cells have been made using Cu(InGa)Se2 films grown by co-deposition from elemental sources.2 The approach considered in the present work is the manufacture of Cu(InGa)Se2 films by deposition onto a flexible substrate using a roll-toroll processing scheme. The film is deposited by thermal evaporation from a series of elemental sources located sequentially through the deposition zone. The process is semicontinuous with interruptions to change rolls. 1.1. Process Description. The process, a schematic of which is shown in Figure 1, can be split into two subsystems: (1) The effusion subsystem consists of the elemental sources for copper, gallium, and indium, where electric energy to resistive heaters is supplied as the input and the outputs are temperature and absorbance, measured by thermocouples and atomic absorption spectroscopy for each source, and selenium is provided in excess and assumed to react stoichiometrically. (2) The deposition subsystem consists of the * To whom correspondence should be addressed. Tel.: (805)893-8133. Fax: (805)893-4731. E-mail: doyle@engineering. ucsb.edu.
Figure 1. Geometry for deposition model.
moving substrate to which the elemental vapor fluxes are delivered as the input and the outputs are the final film thickness and composition measured by X-ray fluorescence (XRF). The XRF sensor is located in the middle above the substrate where it exits the deposition zone, i.e., on the right-hand side in Figure 1. It was developed at ITN Energy Systems in Littleton, CO, and Global Solar Energy in Tucson, AZ;3 the experimental deposition system described in this paper does not have an XRF sensor installed. 1.2. Process Chemistry. The chemical phases of interest for Cu(InGa)Se2 solar cell production are depicted in the phase diagram in Figure 2 and lie on the line connecting In2Se3 to Cu2Se.4 Solar cell efficiency can be substantially improved by partial substitution of indium with gallium in the InSe sublattice.2 In the stoichiometric formula, this is indicated by GaxIn1-x, where x is the molar gallium ratio
GGI ≡
NGa NGa + NIn
10.1021/ie030136g CCC: $27.50 © 2004 American Chemical Society Published on Web 12/17/2003
(1)
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 567
Figure 2. Schematic of the ternary Cu-In-Se phase diagram.4 The phases of interest lie on the line connecting In2Se3 to Cu2Se. R is the desired Cu(InGa)Se2 phase.
which is a measure of the gallium content. N is the number of moles. The other characteristic quantity is the molar copper ratio
CGI ≡
NCu NGa + NIn
(2) 2. Model
which is a measure of the copper content. It determines which of the phases along the In2Se3-Cu2Se tie-line is formed. For copper ratios less than 1, (GaxIn1-x)2Se3, Cu(GaxIn1-x)5Se8 (γ), Cu(GaxIn1-x)3Se5 (δ), or Cu(In1-xGax)Se2 (R) phases occur. For copper ratios greater than 1, a Cu2Se phase is formed. This copperrich phase is highly conductive and will lead to short circuiting of the solar cell’s front and back contacts (shunting); this renders the solar cell useless. For thin-film solar cell production, R is the desired phase; it is formed for copper ratios from approximately 0.9 to 1. Because the Cu2Se phase is directly adjacent to it, there is a hard constraint in keeping the copper ratio less than 1. The process is designed to deliver controlled amounts of copper, gallium, and indium to the substrate. Selenium vapor, on the other hand, is provided in excess to the reaction environment and is built into the film as required by the stoichiometry. If the selenium ratio is defined as
SCGI ≡
NSe NCu + NGa + NIn
(3)
then its dependence on the copper ratio is given by
SCGI )
1 CGI + 3 2 CGI + 1
the inner loop. Therefore, upon inclusion of the inner control loop, the overall control performance might degrade slightly. Because of the transport delay caused by the moving web, effective control requires time delay compensation. A straightforward approach is implementation of a Smith predictor or an internal model controller (IMC). As the most comprehensive framework for control, model predictive control (MPC) is compared to IMC in this work. It is important to note that the design of IMC or MPC controllers requires a process model. 1.4. Overview. The paper is organized in three main sections: Development of the model comprises four parts: (i) the deposition model, (ii) the XRF sensor model, (iii) the Wiener model structure and decoupling strategy, and (iv) the model uncertainties. The control design section comprises the linear model predictive control framework and the internal model control structure. Finally, case studies of the various controllers are presented and discussed.
(4)
1.3. Control Objective and Structure. A critical issue in commercializing this approach is to control the composition and thickness of the Cu(InGa)Se2 film over long deposition times. Because of long startup and shutdown times, it is equally important to efficiently control set-point changes such that multiple recipes can be tested in a single experiment. Analysis of various controllers for reference tracking performance is therefore the focus of this paper. Given the compartmental structure of the combined evaporation and deposition process, a cascaded control structure is proposed. The outer control loop regulates the deposition process by providing flow-rate set-points for the inner effusion control loop. The objective of this paper is to find a suitable control algorithm for the outer loop. It assumes perfect and instantaneous control of
A fundamental model for this process is developed to better understand the physical system and to predict the film growth. The latter is important for controlling the thickness and composition of the deposited film. Modeling of the deposition subsystem requires prediction of the spatial distribution of material on the moving substrate as a function of time. The XRF sensor is modeled as a nonlinear map from the numbers of moles of the individual elements to the characteristic data of the final film: CGI, GGI, and thickness. The model can be cast in a Wiener model structure5 consisting of a linear state space formulation describing the deposition of material onto the substrate, followed by a nonlinear static map describing the XRF sensor. The model development is completed by an uncertainty description. Although most uncertainties are directly related to the deposition subsystem, some address the connection to the effusion subsystem and to its delivery of material to the moving substrate. This section is split into four parts: development of (i) the deposition model and (ii) the XRF sensor model, (iii) introduction of the Wiener model structure and decoupling strategy, and (iv) discussion of model uncertainties. 2.1. Deposition Model. 2.1.1. Spatial Flux Distribution. The local flux on the web surface depends on the spatial flux distribution of vapor leaving the nozzle. To properly describe the flux distribution, a coordinate system is required. Whereas the deposition of material onto the web is most easily described by a Cartesian coordinate system, the flux distribution from each individual nozzle is most conveniently captured by a spherical coordinate system with its origin located at the center of the nozzle exit (Figure 3, left). The coordinates of point P are the zenith θ, the azimuth φ, and the radius r. Because the flux profile is symmetric about the z axis, the azimuth is not required. The geometry of the flux distribution is shown on the right side of Figure 3. The elliptical beam shape has its maximum right above the nozzle and is zero in the nozzle plane. The local flux at the surface element dAn (normal to the beam) depends on the zenith, the radius, and the total flow m ˘ vap leaving the nozzle. The theoretical result for an isothermal cell with an infinitesimally
568
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
Figure 3. Geometry for local flux distribution: spherical coordinate system (left) and elliptical beam shape (right).
small opening (Knudsen cell) is
dm ˘ 1 cos(θ) ) m ˘ vap dAn π r 2
(5)
which is known as the “cosine law of emission”.6 The area element of interest is dA; its orientation relative to the nozzle plane is R. Substituting dAn ) dA cos R into eq 5, the final formulation for the Knudsen cell becomes
dm ˘ 1 cos(θ) cos R m ˘ vap ) dA π r2
(6)
A more realistic beam distribution is described by a cosn(θ) rather than by a cos(θ) functionality.7 Here, n is the beam collimation parameter; the larger n, the more focused the beam. Physically, n is related to the nozzle aspect ratio Λ ) Lnozz/Rnozz, with larger n for larger Λ. The material deposited per unit area then becomes8 n
dm ˘ n + 1 cos (θ) cos R m ˘ vap ) dA 2π r2
(7)
For the sources employed in the experimental facility, a beam collimation parameter of n ≈ 3 was determined experimentally by placing a semicircular hoop over one nozzle and then depositing material onto small pieces of substrate arranged on the hoops perimeter (see Jackson et al.9 for more details on the experimental setup). 2.1.2. Simplifying Assumptions. A rigorous model of the deposition process would describe the local compositions of copper, gallium, indium, and selenium, as well as their reaction products, along, across, and throughout the deposited film as a function of time. In addition to the actual deposition and re-evaporation of material, reaction kinetics and diffusion of material would have to be considered. Because the model is primarily intended for on-line control, several simplifications are made: 1. The growth processes of copper, gallium, and indium are modeled independently, i.e., the reaction kinetics are neglected. This is possible because the XRF sensor is able to measure these elements independently. 2. The layering of material throughout the thickness of the film and the interdiffusion of material is neglected for the same reasons. 3. Diffusion along and across the film is neglected because there are no strong gradients in these directions. 4. The fraction of material that sticks to the substrate (sticking coefficient6) is assumed to be unity. The true
value of this parameter has to be determined experimentally; however, it is of no importance for the current study. Uncertainty in the sticking coefficient is covered by considering the uncertainty in the total flow rate of material delivered to the substrate. 5. Selenium, which is provided in excess, is assumed to be built into the film stoichiometrically on the basis of the phases formed between the two limiting cases (GaxIn1-x)2Se3 and Cu2Se. The stoichiometric amount of selenium is given by eq 4. An interesting extension of the present work, suitable for detailed design studies, could be obtained by removing the first two assumptions to predict a depth profile of material throughout the film. 2.1.3. Prediction of Local Flux. Using the simplifying assumptions discussed above, the crucial step in predicting the film growth is to determine the local flux at any point across or along the moving web as a function of time. For a symmetric setup of the effusion sources under the web, only one-half of the profile (e.g., from the web center to the front) is required; however, it is desirable to keep the prediction as general as possible (e.g., to simulate the effects of asymmetric nozzle placement). The origin of the global Cartesian coordinate system is located on the front left corner of the nozzle plane (see Figure 1). The web has a width of W, a length of L, and a distance from the nozzle plane of H, and it moves at a constant speed v in the positive x direction. The web is discretized into Nx finite elements of size ∆x along its length (machine direction) and into Ny finite elements of size ∆y across its width (cross direction)
∆x )
L Nx
∆y )
W Ny
(8)
Whereas the coordinate system is global to the deposition system, the individual flux distributions are defined locally to each nozzle as discussed in section 2.1.1. The easiest way to compute the local flux at any coordinate (x, y) is to define a local (ξ, ζ) coordinate system for every nozzle (xnozz, ynozz) as shown in Figure 1
ξ ) x - xnozz ζ ) y - ynozz
(9)
For the current geometry, the local flux as given by eq 7 can be simplified, given that the web plane is parallel to the nozzle plane, i.e., R is equal to θ n+1
dm ˘ n + 1 cos (θ) ) m ˘ vap dA 2π r2
(10)
Now, the local spherical coordinates must be converted into local Cartesian coordinates. The radii r and F follow from the Pythagorean theorem as
F ) xξ 2 + ζ 2
(11)
r ) xH 2 + F2 ) xH 2 + ξ 2 + ζ 2
(12)
The value of cos(θ) follows from its definition as
[
]
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 569
cos(θ) )
H H ) r 2 xH + ξ 2 + ζ 2
(13) M≡
The final formulation for the local flux then becomes
n + 1 H n+1 dm ˘ ) m ˘ dA 2π r n+3 vap
(14)
For n ) 3, this simplifies to
H4 2 dm ˘ ) m ˘ dA π (H 2 + ξ 2 + ζ 2)3 vap
(15)
2.1.4. Cinematic Model Formulation. For the solution of the mass balance equations on the moving web, two different numerical techniques were employed. In the finite difference approach, the film is discretized spatially while a set of coupled ordinary differential equations is solved in the time domain. In the cinematic technique,10 discretization in both space and time is performed. The main advantage of the finite difference approach is that it is a continuous technique and, therefore, is better suited for including mass transfer or reaction kinetics. A disadvantage is the computational cost of solving many coupled ODEs. The cinematic approach, on the other hand, is a discrete technique without coupling. A potential disadvantage is that on-line changes of web speed cannot be performed without changing either the sampling time or the number of finite elements. Because changes of web speed are not desired and neither diffusion nor reaction are included in the model, the cinematic technique is preferred for its higher computational efficiency. In the cinematic approach, the movement of the web is simulated step by step. Computing the mass balance reduces to accounting for the newly deposited mass to every finite element during every time step. The web is then moved one ∆x step in the positive x direction. The first row enters on the left with zero mass, and the last row exits on the right. The time step size depends on the web speed and the discretization in the x direction according to
∆t )
∆x v
(16)
The amount of mass ∆m deposited in every time step to every finite element follows from eq 14. Because the actual mass deposited and not the local flux during that time step is needed, this equation has to be multiplied by ∆t and the size of the finite area element ∆A
∆m )
n + 1 Hn+1 m ˘ ∆A∆t 2π rn+3 vap
(17)
To facilitate the algorithm, a shifting mechanism is needed. Starting with all finite elements initially being empty, at every subsequent time step, an “empty” row enters on the left and a “full” row exits on the right. A straightforward way to do this is to organize the finite elements in a mass matrix M according to their (x, y) indices. A web with Nx elements in the machine direction and Ny elements in the cross direction is represented by
m1,1
· · ·
m1,Ny
· · ·
··
· · ·
·
mNx,1 · · · mNx,Ny
(18)
such that the finite elements in the machine direction are the columns and the finite elements in the cross direction are the rows of M. The shifting should be such that the new first row contains 0’s and the last row is either discarded or recorded for future reference. It is accomplished by multiplication by an appropriately chosen shifting matrix S, which is a square matrix of dimension Nx with 1’s on the diagonal below the main diagonal. Here, for example, is a 3 × 2 case
where k indicates the sampling index. With the shifting mechanism, prediction of the film growth becomes straightforward. Initially, all finite elements are zero. The simulation consists of six steps: (i) perform the shifting operation; (ii) loop over all nozzles, all x, and all y to compute the local mass gain; (iii) add the local mass gain to the appropriate finite element; (iv) compute the Se content according to the Cu, Ga, and In contents for the exit row; (v) compute CGI, GGI, and the film thickness for the exit; and (vi) compute the number of moles per unit area, N ′′, at the exit and output the results, and then repeat from step i. In a more compact mathematical notation, the algorithm is written as
Mk ) shift(Mk-1) + ∆mk-1 N ′′k )
1 Mexit,center M ∆A k
(20)
where M is the molecular weight. Although this algorithm works fine, it is computationally inefficient because of the loops over all nozzles, all x, and all y coordinates when computing the mass gain. The next section describes a more efficient way to perform this computation. 2.1.5. State Space Formulation. Having developed the basic algorithm, the cinematic model will now be cast into a discrete state space formulation. Not only does this facilitate the application of standard design and analysis tools for linear systems, it also significantly increases the computational speed. The discrete state space form is given by
xk ) Axk-1 + Buk-1 yk ) Cxk + Duk
(21)
By comparison to the terms in eq 20, the state vector is given by the mass matrix, the shifting operation S forms the A matrix, the vapor flow m ˘ vap forms the input, the local mass gain is described by the B matrix, the
570
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
computation of the number of moles per unit area at the exit is described by the C matrix, and the D matrix is zero. The specific details are summarized below. State Vector. The state vector is directly related to the state matrix defined in eq 18. The most straightforward way of defining the state vector is to stack the state matrix elements columnwise such that
sought-for element has the indices [Nx, Ny /2]. For odd Ny, round up to the next full index
x ) (m1,1, ..., mNx,1|‚‚‚|m1,Ny, ..., mNx,Ny)T
The C matrix, similarly to the A matrix, is very sparse. Its dimension is [1 × (NxNy)]. D Matrix. Although the D matrix is zero, its dimension has to match the dimensions of the other matrices. Because u is [2 × 1] and y is [1 × 1], D is given by
(22)
where, for better legibility, vertical bars have been added to separate the (original) columns. The dimension of the state vector is [(NxNy) × 1]. A Matrix. Because the state matrix has been stacked by columns, the A matrix becomes a block diagonal matrix where the nonzero blocks are formed by the shifting matrix S. For illustration, take the 3 × 2 example from eq 19, which transforms to
[ ]
S Z x Z S k-1
xk )
(23)
where Z is a zero matrix of appropriate dimension. The dimension of the A matrix is [(NxNy) × (NxNy)]. From a computational perspective, it is important to realize that the A matrix is very sparse. Its density, i.e., the number of nonzero elements divided by the number of all elements, is only 0.98%. When implementing the state-space-based algorithm, performance can be increased by exploiting this sparsity. Input Vector. The mass on the web changes through the addition of mass at the flow rate m ˘ vap. Because each source has two nozzles, the input vector for each element is given by
u)
[][ ] u
F
F
)
B
u
m ˘ vap m ˘ Bvap
(24)
where F and B refer to the front and back nozzles, respectively. The dimension of the u vector is [2 × 1]. B Matrix. Each finite element, i.e., each state, is influenced differently by the these inputs. This local gain is described by the B matrix, whose elements, by way of comparison with eq 17, are given by
b
F,B
n+1
n+1H ) ∆A∆t 2π r n+3
(25)
To form the B matrix, these gain matrix elements, as for the state matrix, are stacked columnwise into two columns
)
(
B ) (BF, BB)
F F F F b1,1 , ..., bN |‚‚‚|b1,N , ..., bN x,1 y x,Ny B B B B b1,1 , ..., bN |‚‚‚|b1,N , ..., bN x,1 y x,Ny
)
T
(26)
where the dimension of the BF and BB matrices is [(NxNy) × 1] and the dimension of the B matrix is [(NxNy) × 2]. C Matrix. The output of the deposition model is the scalar number of moles per unit area of the exit center finite element. Thus, the C matrix has to pick out that element and convert from mass to number of moles per unit area via division by M∆A. In the state matrix, the
D ) [0, 0]
(28)
2.1.6. Overall Deposition Model. The results thus far consider one element at a time; for instance, the copper submodel describes the effect of the two copper flow rates on the number of moles of copper per unit area at the film exit. The overall deposition system, on the other hand, consists of copper, gallium, indium, and selenium. To complete the model, the three state space models for copper, gallium, and indium need to be combined to form the full deposition model. Finally, the output has to be augmented by the selenium content, which, according to eq 4, can be rewritten as
y Se )
yCu + 3( yGa + y In) 2
(29)
The complete deposition model then becomes
()[ xCu xGa x In
()[ yCu yGa y In y Se
]( ) [ ]( ) ]( )
0 ACu 0 xCu ) 0 AGa 0 ‚ xGa 0 0 AIn x In k
CCu 0 ) 0 1 Cu C k 2
k-1
0 BCu 0 uCu Ga + 0 B 0 ‚ uGa 0 0 BIn u In 0
k-1
0 0
xCu C Ga 0 C In ‚ x 3 Ga 3 In x In C C 2 2 Ga
(30)
k
2.2. XRF Sensor Model. The deposition system is equipped with an XRF sensor that measures the numbers of moles of Cu, Ga, In, and Se. For better understanding of the sensor uncertainties, the principle of operation is briefly reviewed.3 The sensor consists of two parts, an emitter and a receiver. The emitter sends an X-ray beam onto a designated location of the film, and the receiver measures the energy and count rate of the fluoresced X-rays. The measurement noise is inversely proportional to the square root of the counts. Because of fluorescence from the thicker substrate, there is a tradeoff in signal-to-noise ratio on one side and sampling time and X-ray tube current on the other. Experimentally, a sampling time of 2 min seems to give the best results. A complicated model is then employed to relate the X-ray counts to the individual compositions. Finally,
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 571
CGI, GGI, and thickness are computed and provided as a measurement. The model used here is simplified by the fact that the X-ray fluorescence need not be modeled. The compositions follow directly from the mass per finite element. The model consists of two parts: (i) averaging the measurement over the sample period of 2 min and (ii) converting from composition, i.e., number of moles per unit area, to CGI, GGI, and thickness. The first part is most easily implemented by integrating the process model output over 2 min, holding that value for 2 min, and finally dividing by 2 min to obtain the average. The integrator is then reset, and the measurement cycle starts over again. The second part results in a straightforward nonlinear map. In addition to eqs 1-3, an equation describing the thickness H is needed
H)
mCIGS
)
∆AFCIGS
∑ im i
)
∑ iN ′′i Mi FCIGS
∆AFCIGS
(31)
where N ′′ represents the number of moles per unit area, the output from the linear deposition model; FCIGS is the overall film density; and i ) Cu, Ga, In, Se. For ease of notation, the following abbreviations are introduced: Following the standard notation, the outputs from the linear deposition model are called y, i.e., (y1, y2, y3, y4) ≡ (N ′′Cu, N ′′Ga, N ′′In, N ′′Se). The outputs from the XRF sensor are called z, i.e., (z1, z2, z3, z4) ≡ (CGI, GGI, H, SCGI). The next step involves an inversion such that yi can be computed as a function of zi. This inverse map will be needed in section 3 for decoupling the control loops. Using these abbreviations, the full forward and reverse transformations are given by
y1 ) y2 )
y3 )
y4 )
z1z3
z1 )
D z2z3
z2 )
D
(1 + z1)z3z4 D
y2 + y3 y2 y2 + y3
∑i)1yiMi 4
(1 - z2)z3 D
y1
z3 ) z4 )
(32)
FCIGS y4 y1 + y2 + y3
where D is the common denominator of all inverted equations
D)
z1 M1 + z2 M2 + (1 - z2)M3 + (1 + z1)z4 M4 FCIGS
(33)
2.3. Wiener Model. The deposition and XRF models can be combined to form a Wiener model.5 The unique characteristic of a Wiener model is that of a linear dynamic element followed by a static nonlinear element (see Figure 4). For process control (see section 3), invertibility of the XRF model can be used to remove the nonlinearity from the control problem. In the current case, this proves particularly useful, because the linear model is de-
Figure 4. Wiener model structure: a linear dynamic element, here the state space deposition model, followed by a static nonlinear element, here the XRF model.
coupled, i.e., each flow rate directly influences the corresponding number of moles per unit area at the deposition exit. 2.4. Uncertainty Description. The model is completed by a discussion of the uncertainties, which can be categorized as input uncertainties, model uncertainties, and output uncertainties. They are caused by unmodeled dynamics, simplifying model assumptions, parametric uncertainty, and noise. 2.4.1. Input Uncertainty. There are two kinds of input uncertainty, total flow and flow distribution. Both are related to the effusion subsystem because of its delivery of material to the moving substrate. Total Flow. The first contribution is uncertainty in the total flow due to the fact that the flow will be predicted by the effusion model. Thus, uncertainties in the flow prediction can result in a deviation of the actual flow from the predicted flow. In the simulation studies, this has been taken into account by a steady offset of the actual flow rate. Another effect, which has not been analyzed here, is a monotonic decay in flow rate due to nozzle clogging. Flow Distribution. The second contribution is uncertainty in the distribution of the flow between the front and back nozzles. In the ideal case, the total flow from the effusion cell is split evenly between the two nozzles; this is based on the assumption that the melt has a uniform temperature. If gradients were present, for instance, the front hotter than the back, the total flow would be split unevenly, for instance, 60% in the front and 40% in the back. For the experimental system considered in this study, which has the nozzles arranged symmetrically under the web, this uncertainty has no effect on the output, as long as the measurement is taken in the middle (symmetry) and assuming that both the front and back nozzles have the same beam collimation parameter. 2.4.2. Deposition Model Uncertainty. The deposition model’s uncertainties are both of parametric nature: beam collimation and time delay. Beam Collimation. The beam collimation parameter was determined by fitting the beam distribution eq 7 to experimental data. The value of n0 ) 3 has an uncertainty of approximately 10%. When testing the controller, this uncertainty can be taken into account by using the nominal value n0 in the controller model and, for instance, a value of n ) 3.3 in the process model. A more realistic uncertainty analysis would consider a dynamically changing value of n, for instance, due to nozzle clogging. Then, its effect on the B matrix in eq 21 needs to be quantified
xk ) Axk-1 + B(nF,nB)uk-1
(34)
Linearization in a Taylor series about the operating point OP ) (x0, u0; nF0, nB0) yields
(
n -n xk ) Axk-1 + Buk-1 + [BFn uF0 , BBn uB0 ] nF - nF0 B B0
)
(35)
where B ) B(nF0,nB0), BFn ) BFn (nF0), and BBn ) BBn (nB0). For brevity of notation, the arguments are neglected
572
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
with the understanding that both B and Bn are always computed for their nominal parameter values. The elements bF,B of B are given by eq 25, and the elements of Bn ≡ ∂B/∂n are
bF,B n )
1 H ∂bF,B ) bF,B + ln ∆A∆t ∂n n+1 r
[
( )]
(36)
Time Delay. The cinematic model is a simple gain (B matrix) plus time delay (A matrix) model; the higher the discretization Nx in the machine direction, the more precise the modeling of the time delay. This results in different types of uncertainty: If the process model is not computed for a fine discretization (Nx > 50), the simulation is not an accurate representation of the experimental system’s dynamics. If a state space model is used for controller design, a low-order model results in sluggish control. On the other hand, a controller based on a high-order model yields theoretically perfect control. If a first-order plus time delay model is used for controller design, the ideal pure delay model is approximated by a time constant and regular time delay, e.g., G(s) ) [K exp(-Rs)]/(Ts + 1). 2.4.3 Output Uncertainty. The output uncertainty is due to measurement noise and parametric uncertainty in the XRF sensor model. Measurement Noise. As discussed in section 2.2, the XRF measurement noise depends on the count rate and therefore on the sampling time. The research group at ITN in Littleton, CO, determined typical levels of measurement noise for a sampling time of 2 min. These values are reported in Table 1. Table 1. XRF Noise element
noise (%)
Cu Ga In Se
0.68 0.89 1.7 0.5
Composite Film Density. As discussed in section 2.2, the XRF sensor model contains the composite film density FCIGS as a parameter. The nominal value for pure Cu(InGa)Se2 is approximately 5.7 g/cm3. Depending on the composition, the film can contain any two adjacent of the five phases shown in Figure 2. The density of the desired Cu(In1-xGax)Se2 phase depends on the gallium ratio x. For the two limiting cases x ) 1 (CuGaSe2) and x ) 0 (CuInSe2), values of 5.45 and 5.65 g/cm3, respectively, are reported in the literature.11 This is a variation of less than 5% in the nominal value. Although data for the other phases are not available in the literature, an uncertainty of no more than (10% seems appropriate. 3. Control Design The inputs and outputs of the combined linear deposition and nonlinear XRF models are highly coupled with two of the measured outputs (CGI, thickness) depending on all three vapor fluxes, and the third output (GGI) depending on two vapor fluxes. Controller performance is greatly improved by decoupling the individual control loops through inversion of the nonlinear map and controlling the linear deposition system
Figure 5. MPC structure with feedback via bias term. The simulation model is given by its Wiener model structure LIN f NL. Decoupling is achieved via NL-1, i.e., by inversion of the nonlinear map. The noisy output is filtered by DF. MPC output correction is achieved via the bias term yB.
directly. Because of the nonlinear model structure of the combined model, this is easily accomplished. The decoupled system can then be controlled by three independent control loops for the individual effusion sources. An apparent alternative to full decoupling is control of the film thickness by manipulation of the web speed. This approach is inherently complicated given that the cinematic model depends on the web speed as described by eq 16. Because the individual effusion rates need to be manipulated regardless, direct control of the three sources via decoupling is a more elegant solution. Because of the transport delay caused by the moving web, effective control requires time delay compensation. A straightforward method of accomplishing this is implementation of a Smith predictor in combination with a PI controller. However, because the goal of this paper is the comparative study of different advanced controllers, this is not our first choice. The most comprehensive framework for the control of constrained multivariable processes is model predictive control (MPC). The Wiener model structure can be used to decouple and linearize the control loops12 such that the well-established linear MPC technique13 is sufficient for control of the deposition system. Because the decoupled deposition process is only single-input-single-output with simple constraints (flow rates cannot be negative), internal model control can potentially combine the advantages of unconstrained MPC, but avoid its disadvantages. For instance, MPC allows for easy on-line tuning via adjustment of physically meaningful parameters without any concern about closed-loop stability and for good performance without intersample rippling.13 3.1. Model Predictive Control. The MPC algorithm used in this paper combines the Wiener-model-based decoupling of Norquay et al.12 with the state-spacemodel-based output prediction of Lee et al.14 Regarding feedback, two choices can be made: (i) the difference between model prediction and measurement can be used as a bias to predict future outputs,13 or (ii) controller performance and disturbance rejection can be improved by state estimation, for instance, by means of a Kalman estimator14 and augmentation of the model by a disturbance model.15 3.1.1. MPC with Bias Term. The basic MPC algorithm uses the difference between model prediction and measurement as a constant bias to correct the prediction of future outputs.13 The control flowsheet (Figure 5) has been augmented by a digital filter block DF to reduce the XRF measurements noise. A Butterworth filter was selected as it provides the best Taylor series approximation to the ideal low-pass filter response at analog frequencies and the magnitude response is maximally flat in the passband and monotonic overall.16
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 573
Figure 6. MPC structure with state estimation. The simulation model is given by its Wiener model structure LIN f NL. Decoupling is achieved via NL-1, i.e., by inversion of the nonlinear map. The states are corrected via the Kalman filter KF.
This approach is limited by two factors: First, it is unlikely that the bias stays constant over the prediction horizon, which might result in either too aggressive or too sluggish a control action. Second, the controlled outputs need to be measurable. Both of these disadvantages can be avoided by state-estimation-based MPC. 3.1.2. MPC with State Estimation. In this approach, state correction is accomplished by means of a Kalman estimator. Although the correction occurs only at the current sampling time, it is not constant over the prediction horizon because x subsequently is multiplied by A, A2, and so forth. This setup (Figure 6) does not require a digital noise filter as the Kalman filter not only estimates the disturbance, but also filters the measurement noise. Correction of the states can be greatly improved by augmentation of the original model with a disturbance model. Disturbance Model. Process knowledge is useful in devising a good disturbance model; this information will be reflected in the disturbance matrix Bd. It is common to model the disturbance as integrated white noise via identity state matrices Aw, Bw, and Cw.17 The general form of the augmented model is then
xk ) Axk-1 + Buk-1 + Bddk-1
(37)
yk ) Cxk + vk-1
(38)
w w w xw k ) A xk-1 + B wk-1
(39)
dk ) Cwxw k + wk-1
(40)
where wk and vk are white noise. In the given case, the effects of disturbances in the input and the beam collimation on the model states are known, namely, via the B and the Bn matrices. Under a [2 × 1] input disturbance du, eq 37 takes the following form u xk ) Axk-1 + B(uk-1 + d k-1 )
(41)
The effect of a [2 × 1] beam collimation disturbance dn is given by eq 35, such that eq 37 becomes n xk ) Axk-1 + Buk-1 + [BFn uF0 , BBn uB0 ]d k-1
(42)
Because both du and dn are unknown disturbances of the same structure, they can be lumped together. Matrix Bd in eq 37 is then given by
Bd ) B + [BFn uF0 , BBn uB0 ]
(43)
Because of the symmetry of both B and [BFn uF0 , BBn uB0 ], there is another degree of redundancy. It therefore suffices to use only, for instance, the B submatrices of the front nozzle and a single [1 × 1] disturbance
Bd ) BF + BFn uF0
(44)
Regarding the output disturbances, the XRF measurement noise is captured by noise term vk, and the density uncertainty is not modeled here because it is not part of the linear state space model. The disturbance model thus far has only a single disturbance state such that the matrices Aw, Bw, and Cw are all unity. The augmented system can be cast into the following form
x k′ ) Φx k-1 ′ + Buk-1 + Γ wwk-1
(45)
yk ) Ξ x′k + vk-1 where
(
)
x ′ ) (x, xw),
Bd A Φ) Z 1,Nx 1
Ξ ) (C, 0),
ZN ,N ZNx,1 Γw ) Z x x 1 1,Nx
(
)
where Za,b is an [a × b] zero matrix. Observability. A prerequisite for the design of a Kalman estimator is that (i) the states of the original model be at least partially observable and (ii) the disturbance model states be fully observable. The condition for observability is that the observability matrix has full rank
Ob ) (C,CA, ..., CAnx-1)T
(46)
where nx is the number of states. If R ) rank(Ob) < nx, then the number of unobservable states is nx - R. First, the observability of the original state space system is checked. Independently of the choices of Nx and Ny, exactly Nx states are always observable. Via the observability staircase form,18 it can be proven that these are the Nx states corresponding to the finite elements in the center row where the XRF sensor is located. If states are not observable, it must be verified that they are detectable, to ensure that they do not grow unobserved in an unstable fashion.19 Alternatively, and this is the approach chosen here, the unobservable states can be removed from the problem by use of a minimal realization. The state space model then only contains the Nx observable states, which has the added benefit of increased computational speed. This minimal realization will be augmented with the disturbance model as discussed above. It turns out to be fully observable, such that the Kalman estimator can be designed next. Kalman Estimator. The optimal Kalman filter gain K is computed via the MATLAB “dlqe” command. The filter is tuned via the state covariance matrix while the measurement covariance is fixed by the XRF sensor noise. The state covariance of the original model is fixed at 10-14, as all model uncertainty is to be covered by the disturbance model. In the tuning of the Kalman filter, there is a tradeoff between disturbance estimation and measurement noise filtering: the covariance of the disturbance states of all three control loops was selected at one order of magnitude below the measurement covariance, thus emphasizing noise filtering. 3.1.3 Controller Tuning. Tuning of the MPC controllers requires selection of the prediction and move horizons, as well as of the input penalty.
574
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
Figure 7. Decoupled IMC structure. The simulation model is given by its Wiener model structure LIN f NL. Decoupling is achieved via NL-1, i.e., by inversion of the nonlinear map. The XRF sensor sampling is modeled by DM, and the noisy output is filtered by DF.
The prediction horizon has to be long enough for effects of an input change on the output to become steady. For the deposition system, the prediction horizon should, therefore, be at least equal to Nx, which models the time delay. For the case studies presented in section 4, Nx ) 50 was used, such that p ) 50 as well. Because the deposition system model contains no dynamics other than the time delay, a move horizon of m ) 1 suffices. The emphasis of this paper is to study the fastest possible reference tracking. Therefore, the input penalty was set to 0. This makes the controller very aggressive; however, it does not become oscillatory or unstable. 3.2. Internal Model Control. 3.2.1. Decoupled IMC Structure. A standard technique for realizing the IMC controller in a classic feedback structure involves a rerouting of the model output to the input of the control block.20 The IMC structure (see Figure 7) can be decoupled through inversion of the nonlinear XRF map, i.e., by using the same decoupling strategy as for MPC.12 Two more elements are part of the control loop: As discussed in section 2.2, the XRF measurement is taken every 2 min as an averaged integral. Given that the linear dynamic model runs at a faster sampling time because of its higher discretization, a digital measurement block DM has been included. Similarly to MPC with bias term (see section 3.1.1), a digital Butterworth filter is used to reduce XRF measurements noise. The best possible internal model controller for this process is based on the full discrete state space model. However, because it is unlikely that such a precise model is available in practice, it is useful to identify simpler models that could also be obtained experimentally on a physical system. A first-order plus time delay model was obtained via step tests on the process model. The resulting control structure is equivalent to a Smith predictor in combination with a PI controller. The parameters of this simplified model were obtained by nonlinear optimization based on the sum-squared error between the step response and prediction. 3.2.2. Controller Design and Tuning. IMC controller design and tuning consists of two parts. For design, the process model has to be inverted where possible.20 The tuning is accomplished by way of a filter.21,22 The filter equally serves to render the controller at least semiproper.23 If the inverted model is already semiproper, a filter should still be included for tuning. State Space Controller Model. The discrete state space model is a pure gain (B matrix) plus time delay (A matrix) model. Because the time delay cannot be inverted, the controller consists of the inverse gain KSS plus a first-order digital IMC filter
CSS )
1 λz - 1 KSS λz z - 1
(47)
where λz > 1 is the filter tuning parameter. The closer λz is to 1, the more sluggish the controller is. Here, a value of λz ) 3 was selected. This makes the controller very aggressive; however, it does not become oscillatory or unstable. First-Order Plus Time Delay Controller Model. The continuous model is fully invertible except for the time delay. A first-order filter is needed to make the controller semiproper
Cstep )
1 Ts + 1 Kstep λss + 1
(48)
The filter tuning parameter λs was chosen identical to the model’s time constant. This achieved the best performance with minimal oscillation. Smaller values of λs result in very oscillatory behavior; larger values in sluggish response. 4. Case Studies Simulation results are presented for the four control structures discussed. For objective comparison with MPC, the digital measurement block has been removed from the IMC structures. For accurate representation of the time delay, both process and controller models were discretized with Nx ) 50 finite elements along the machine direction. Controller performance is compared for reference tracking under various input disturbances, parametric uncertainty, and measurement noise. For Cu(InGa)Se2-based solar cells, a typical operating point of Z0 ) (0.85, 0.4, 2 µm) was selected. The operating point in Y is then Y0 ) (3.20, 1.51, 2.26) µmol/ cm2. Finally, this operating point can be achieved by the following steady-state flow rates: U0 ) (0.72, 0.34, 0.92) g/h‚nozzle. As a reference trajectory, a step of +10% from the operating point of CGI and GGI at constant thickness was chosen. For input uncertainty, a steady offset of the flow rate was assumed. The three values [0, -3, -10]% were tested. For beam collimation and film density, an uncertainty in each of +10% of the nominal value was assumed. The results given in Figure 8 are based on noise-free data; the results in Figure 9 are based on noisy data according to the XRF measurement noise given in Table 1. The substrate moves at a speed of v ) 1 in./min through the deposition zone. With the length of the deposition zone of L ) 15 in., the minimum time required for a set-point change is 15 min. This time will henceforth be referred to as a web pass. Time Response Analysis of Noise-Free Data. Figure 8 shows the time response of the copper control for 0, -3 and -10% input disturbances. All controllers are tuned for fast reference tracking. Note that the input uCu for the nominal output yCu,0 before the step is not zero as would be typical for a normalized deviation variable. The cause for this is the parametric uncertainty of beam collimation and film density. At low levels of input disturbance, IMC based on a first-order plus time delay model takes about one web pass longer to reach the new set-point than IMC based on a state space model, whereas the two types of MPC show very similar responses. At high levels of input disturbance, all controllers need about one extra web pass to reach the set-point. Relative to each other, the two IMC-based controllers show the same behavior as at low input disturbance. For MPC, on the other hand,
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004 575
Figure 8. Comparison of copper control time responses for three different input disturbances: input u-Cu is the normalized Cu flow rate, and output n-Cu (NCu ′′ ) is the normalized number of moles per unit area. Legend: IMC with first-order plus time delay model (thin solid), IMC with state space model (thin dotted), MPC with bias term (bold dotted), MPC with state estimation (bold solid), set-point (dash-dotted).
now, the controller with state estimation is about one web pass faster than the controller based on output bias. Mean and Root-Mean-Square Analysis of Noisy Data. To analyze controller performance of noisy data, the simulation data are processed as follows: First, the means and standard deviations of the filtered outputs Y1, Y2, Y3, Z1, Z2, and Z3 are determined at every sampling time k for N ) 100 simulations; for instance, for Y1(i,k), where i refers to the ith simulation
Y1(k) )
σY1(k) )
{
1
1
N
∑ Y1(i, k) N i)1
∑
}
0.5
N
N - 1 i)1
(49)
[Y1(i, k) - Y1(k)]2
(50)
Using these results, three characteristic scalar quantities are obtained by averaging over the simulation time horizon k ) k1 ... k2, where k1 designates the sampling index one web pass after the reference change and k2 the sampling index just before a new reference change. The mean tracking error (Figure 9, row 1)
MEAN(eY1) )
1
k2
∑ [Y1(k) - Y d1]
k2 - k1k)k1
(51)
where Y d1 is the set-point, is a measure for how closely on average the controllers track the set-point and how large the average deviation is from the set-point. With an average deviation of less than (0.4, all four controllers track perfectly for all input disturbances. Also, the results from worst to best are in the order IMC with step response model, IMC with state space model, and MPC. These observations are confirmed by the analysis of the root-mean-square (RMS) deviation of the tracking
Figure 9. Controller comparison for three input disturbances. The simulation data are processed as follows: First, the means and standard deviations of Y1, Y2, Y3, Z1, Z2, and Z3 are determined at every sampling time for 100 simulations with filtered data. Then, a scalar quantity for each variable is obtained by averaging these data over the simulation time horizon (see text for details). Row 1 displays a plot of the mean tracking error; row 2, the rootmean-square (RMS) tracking error; row 3, the RMS standard deviation; and row 4, the same as row 3 but using unfiltered data. Legend: IMC with first-order plus time delay model (white), IMC with state space model (dark gray), MPC with bias term (light gray), MPC with state estimation (black).
error (Figure 9, row 2), which is defined as
RMS(eY1) )
{
}
k2
1
∑
0.5
k2 - k1 k)k1
[Y1(k) - Y d1] 2
(52)
and is a measure for the absolute deviation from the set-point that penalizes large errors. Analysis of the RMS is important because the mean value can be misleading as a result of potential cancelations of positive and negative deviations. Here, it becomes apparent that MPC with state estimation is better than MPC with bias term. Furthermore, whereas the deviation increases with increasing input disturbance, at higher input disturbances, the MPC-based controllers are clearly superior to the IMC-based controllers. Last, the root-mean-square standard deviation (Figure 9, row 3)
RMS(σY1) )
{
1
k2
∑
k2 - k1 k)k1
}
0.5
[σY1(k)]2
(53)
is computed. The filters work effectively for all input disturbances. The noise level does not increase with input disturbance because the noise is caused by the output disturbance. The results for Y1, Y2, and Y3 correspond to the data given in Table 1. Also, to verify the filter performance, this value has been computed and plotted for unfiltered data (Figure 9, row 4). Qualitatively, the results are identical; however, the values are twice as large. Thus, both the Butterworth and Kalman filters effectively reduce the measurement noise. 5. Summary A discrete state space deposition and a nonlinear XRF model for the production of thin-film solar cells via
576
Ind. Eng. Chem. Res., Vol. 43, No. 2, 2004
physical vapor deposition onto a moving substrate have been developed. The combination of these linear and nonlinear submodels comprises a Wiener model. By inverting the nonlinear Wiener model, the structure of the map was utilized in linearizing the problem and decoupling the control loops. Various input, output, and parametric uncertainties were discussed and incorporated in the control studies. Two types of controllers, IMC and MPC, were tested. For IMC, both first-order plus time delay and state space models were used; for MPC, both controllers with output bias and state estimation were developed. Controller performance was compared under input, parametric, and output uncertainty. The MPC-based controllers exhibit shorter rise times and better compensation of input disturbances than the IMC-based controllers. The two IMC controllers perform similarly well, which is remarkable considering that the first-order plus time delay model is much simpler than the state space model. At high input disturbances, MPC based on state estimation does show a significant improvement over MPC with bias term, consistent with intuition. The recommendation for control of a real system thus depends on the level of expected input disturbance. If input disturbances are low, IMC based on an experimentally obtained first-order plus time delay model should be used. On the other hand, if strong input disturbances are expected, MPC with state estimation should be considered. Acknowledgment Partial funding from Global Solar Energy and the Delaware Research Partnership is gratefully acknowledged. Literature Cited (1) Birkmire, R. W. Compound Polycrystalline Solar Cells: Recent Progress and Y2K Perspective. Sol. Energy Mater. Sol. Cells 2001, 65 (1), 17-28. (2) Birkmire, R. W.; Eser, E. Polycrystalline Thin Film Solar Cells: Present Status and Future Potential. Annu. Rev. Mater. Sci. 1997, 27, 625-653. (3) Eisgruber, I. L.; Joshi, B.; Gomez, N.; Britt, J.; Vincent, T. In Situ X-ray Fluorescence Used for Real-time Control of Cu(InxGa1-x)Se2 Thin Film Composition. Thin Solid Films 2002, 408 (1-2), 64-72. (4) Go¨decke, T.; Haalboom, T.; Ernst, F. Phase Equilibria of Cu-In-Se I. Stable States and Nonequilibrium States of the In2Se3-Cu2Se Subsystem. Z. Metallkd. 2000, 91 (8), 622-634.
(5) Pearson, R. K. Discrete-Time Dynamic Models; Oxford University Press: New York, 1999. (6) Maissel, L. I.; Glang, R. Handbook of Thin Film Technology; McGraw-Hill: New York, 1970. (7) Ohring, M. The Material Science of Thin Films; Academic Press: San Diego, CA, 1992. (8) Pulkner, H. K. Coatings on Glass; Elsevier: New York, 1984. (9) Jackson, S. C.; Baron, B. N.; Rocheleau, R. E.; Russell, T. W. F. Molecular Beam Distributions from High Rate Sources. J. Vac. Sci. Technol. A 1985, 3 (5), 1916-1920. (10) Lakshmanan, C. C.; Potter, O. E. Cinematic Modeling of Dynamics of Countercurrent Systems. Comput. Chem. Eng. 1990, 14 (9), 945-956. (11) Hahn, H.; Frank, G.; Klingler, W.; Meyer, A.-D.; Sto¨rger, G. Untersuchungen u¨ber terna¨re Chalkogenide: U ¨ ber einige terna¨re Chalcogenide mit Chalkopyritstruktur. Z. Anorg. Allg. Chem. 1953, 271, 153-170. (12) Norquay, S. J.; Palazoglu, A.; Romagnoli, J. A. Model Predictive Control Based on Wiener Models. Chem. Eng. Sci. 1998, 53 (1), 75-84. (13) Garcia, C. E.; Prett, D. M.; Morari, M. Model Predictive Control: Theory and PracticesA Survey. Automatica 1989, 25 (3), 335-348. (14) Lee, J. H.; Gelormino, M. S.; Morari, M. Model Predictive Control of Multi-rate Sampled Data Systems: A State-space Approach. Int. J. Control 1992, 55 (1), 153-191. (15) Ricker, N. L. Model Predictive Control with State Estimation. Ind. Eng. Chem. Res. 1990, 29 (3), 374-382. (16) MATLAB Signal Processing Toolbox; The MathWorks, Inc.: Natick, MA, 2001. (17) Lee, J. H.; Ricker, N. L. Extended Kalman Filter Based Nonlinear Model Predictive Control. Ind. Eng. Chem. Res. 1994, 33 (6), 1530-1541. (18) Rosenbock, H. H. State-Space and Multivariable Theory; John Wiley & Sons: New York, 1970. (19) Brogan, W. L. Modern Control Theory, 3rd ed.; Prentice Hall: Englewood Cliffs, NJ, 1991. (20) Garcia, C. E.; Morari, M. Internal Model Controls1. Unifying Review and Some New Results. Ind. Eng. Chem. Proc. DD 1982, 21 (2), 308-323. (21) Brosilow, C. B. The Structure and Design of Smith Predictors from the Viewpoint of Inferential Control. In Joint Automatic Control Conference Proceedings; AIChE: New York, 1979. (22) Mehra, R. K.; Rouhani, R.; Rault, A.; Reid, J. G. Model Algorithmic Control: Theoretical Results on Robustness. In Joint Automatic Control Conference Proceedings; AIChE: New York, 1979; pp 387-392. (23) Ogunnaike, B. A.; Ray, W. H. Process Dynamics, Modeling, and Control; Oxford University Press: New York, 1994.
Received for review February 13, 2003 Revised manuscript received July 29, 2003 Accepted October 2, 2003 IE030136G
