Fabrication of Silica Optical Fibers: Optimal Control Problem Solution

In this work, a new approach to solving problems of optimal control of manufacture procedures for the production of silica optical fiber are proposed. The procedure of silica tubes alloying by the Modified Chemical Vapor Deposition (MCVD) method and optical fiber drawing from a preform are considered. The problems of optimal control are presented as problems of controlling distributed systems with objective functionals and controls of different types. Two problems are formulated and solved. The first of them is the problem of the temperature field optimizing in the silica tubes alloying process in controlling the consumption of the oxygen–hydrogen gas mixture (in the one- and two-dimensional statements), the second problem is the geometric optimization of fiber shape in controlling the drawing velocity of the finished fiber. In both problems, while using an analog to the method of Lagrange, the optimality systems in the form of differential problems in partial derivatives are obtained, as well as formulas for finding the optimal control functions in an explicit form. To acquire optimality systems, the qualities of lower semicontinuity, convexity, and objective functional coercivity are applied. The numerical realization of the obtained systems is conducted by using Comsol Multiphysics.


Introduction
Nowadays, fiber optics have solidly taken leading positions among various branches of science and technology [1]. Fiber-optic technologies are widely used in long-distance and local communication systems, instrument engineering, and scientific research [2]. At the same time, plenty of information in this area of science is either insufficiently or fragmentarily studied, mostly out of reach, not reliable enough, or not even reproducible. Such a situation leads us to the fact that to improve the characteristics of the fiber, we have to put out a lot of effort and resources. First of all, this applies to optical fibers, for which there is a need for significant improvement. Thereby, research of the technology features and fiber production is a well-timed and relevant objective.
It is known that the main optical and mechanical characteristics of optical fibers are determined by the structure and quality of the initial billet [3]. However, the determination of these indicators largely depends on the alloying regimes (vapor deposition) and fiber drawing. In the alloying process, such parameters as the pressure and temperature of the reagents at the tube inlet, the speed of the burner, the temperature of the outer surface of the silica tube in the reaction zone, the tube diameter, the tube deflection, and a number of other parameters are controlled. The quality control of preform parameters includes measurements of the refractive index, core diameter, the outer diameter of the preform, and concentricity. In the process of fiber drawing, the speed and strength of the drawing, the outer diameter of the fiber, the concentricity and thickness of the coating, as well as • MCVD is the modified chemical vapor deposition; • OVD is the outside vapor deposition; • VAD is the vapor axial deposition; • plasma-chemical methods (PMCVD, PCVD, etc.).

Motivation
Nowadays, the problems of distributed systems' optimal control gain particular interest because they describe various high-tech production processes, including all the main technological cycles of the silica optical fibers' production. Under the main production cycles, we understand both the process of forming a silica preform (according to MCVD technology for our work) and fiber drawing (the final stage of production). The mathematical modeling of these technological processes is widely researched [10][11][12]. The authors of various literature analyze the features of models that describe the production of different types of special fibers [13,14] and the mutual influence of model parameters; numerical results of these calculations are obtained. The authors of this research did not aim to improve the existing mathematical models or their parameters. The general objective of this work is to restrict the research to the available models that describe the main stages of optical fiber production so we can formulate, justify, obtain, and analyze the results of optimal control problems' numerical solutions for the chosen processes.
The approach proposed in this work is based on two main ideas: (a) To control and manage the allowing processes by the MCVD method and the drawing of optical fibers not by measuring the temperature and fiber diameter at only one point, as it is done now, but by measuring these parameters over extended sections (by length). In other words, to conduct a distributed observation; (b) To control the process not with the help of PID controllers, as it is done now, but based on the theory of the optimal control of distributed systems.

One-Dimensional Mathematical Model of MCVD Process
The main problem is the difficulty in controlling the temperature of the reaction zone, where physical and chemical processes take place. At the same time, the temperature distribution field in the zone of oxide formation is not limited to determining just the size and concentration of the coagulated particles but the direction of their motion under the influence of thermophoresis forces, since the thermophoretic force acting on each particle is proportional to the temperature field gradient at a given point [15]. In MCVD practice, the temperature on the surface of a silica tube is monitored using non-contact infrared pyrometers.
We developed a mathematical model for heating a silica tube with a movable heat source that will describe the temperature field with sufficient accuracy in the silica tube formed during a real technological process. The mathematical model was obtained under the following assumptions: the temperature field of the silica tube is axisymmetric (this is ensured by the rotation of the pipe); the heat exchange with the external environment and gas flowing inside the pipe is described by Newton's law; the radiation from the outer surface of the pipe obeys the Stefan-Boltzmann law.
The model includes the energy equation (in our case, the heat equation) and the model of the moving source, that is, a description of the shape and power of the supplied heat flux, as well as the source motion law [16,17]. The heat flux from the plume changes along the z-axis and is described by the following Gauss function: where v(ξ) is burner moving speed, H is a form parameter (heat, flow, width), q max is the burner intensity (power), t is time, and z is a spatial variable. Then, the heat equation for a silica tube in a cylindrical coordinate system has the following form: where Θ(t, r, z)-temperature; ρ, C p , λ are the density, specific heat, and thermal conductivity of silica, respectively. In general, these characteristics are a function of temperature Θ(t, r, z); r and z are spatial variables, and, moreover r ∈ [r 0 ; R], z ∈ [0; L], and t ∈ [0; T], where r 0 is the pipe's inner radius, R is the pipe's outer radius, L is the length of the pipe, and T is the deposition process duration. The effective thermal conductivity of silica as a material transparent to thermal radiation consists of two components: conductive λ k and radiant λ r . Since the medium is optically thick, the Rosseland approximation was used to estimate the radiation component of the thermal conductivity coefficient as follows: The initial and boundary conditions on the pipe surfaces have the following form: λ ∂Θ(t, r, z) ∂r where α gas and α env are the heat transfer coefficients of gas and the environment (air), respectively, Θ gas and Θ env are gas and ambient temperatures, respectively, ε is the black ratio, and σ 0 is the Stefan-Boltzmann constant. Furthermore, we will omit the arguments when writing functions if this entry is not fundamentally important. Temperature changes along the pipe cross-section are small compared to temperature changes along the z-axis due to the axial symmetry of the temperature field and the small thickness of the pipe wall. That is why we replace the temperature values in the crosssection with the average value over the region D = { (r, φ)|r 0 ≤ r ≤ R; 0 ≤ φ ≤ 2π}. We integrate both sides of Equation (2) over the indicated region as follows: Then, the first integral on the right side of the last equation takes the following form when paying attention to conditions (4) and (5): After averaging the remaining terms over domain D and dividing both sides of the equation by π R 2 − r 2 0 , we obtain the following equation of the following form with properties that depend on temperature: where: (8) is supplemented by the initial and boundary conditions similar to (3), (6), and (7), as follows: If the output of the object can be considered small deviations of ∆Θ(t, z) from some temperature Θ * (t, z) = Θ * , while assuming that Θ(t, z) = Θ * + ∆Θ(t, z), then with fairly smooth dependencies ρ(Θ), C p (Θ), λ(Θ), Equation (8) can be linearized in the vicinity of temperature Θ(t, z)| z=0 = Θ 1 (t), ∂Θ(t,z) ∂z z=L = Θ 2 (t) by expanding nonlinear dependencies into a Taylor series. We use this technique in Equation (8), which is a one-dimensional second-order differential equation with a nonlinear operator calculated as follows: Let this equation also describe the state of the object Θ(t, z) = Θ * + ∆Θ(t, z), where ∆Θ(t, z) is small variation with respect to a certain state Θ * , which corresponds to a certain external influence q * (t, z). At the same time, Θ * and q * (t, z) satisfy Equation (9), that is: Assuming that the function A is twice continuously differentiable by the totality of arguments, we can obtain a linear equation for ∆Θ(t, z) (perform linearization procedure), which is called the equation in variations [18] as follows: These notations are introduced here Θ = ∂Θ ∂t , Θ = ∂ 2 Θ ∂z 2 , Θ = ∂Θ ∂z , where symbol "*" corresponds to the Θ * state of the system.
Derivatives in the Equation (10) present in the following form: Considering the introduced notations, we obtain: Divide both sides of the resulted equation by ρ(Θ * )C p (Θ * ) to obtain: Then, we have the following: where: Note that if the temperature Θ * depends only on the coordinate z, then in expression ∂A ∂Θ * , the first term turns to zero. In the case Θ * = const, then in ∂A ∂Θ * , the first three terms equal to zero and ∂A ∂Θ = 0.. Similarly, we linearize the boundary and initial conditions, as a result of which we obtain the following boundary and initial conditions for the linearized problem: Thus, the linear distributed differential problem (11), (12) is formulated, described by the one-dimensional heat equation, simulating the heat transfer of the MCVD process.

Solution of the Distributed Optimal Control Problem
In general, the goal of influence heat source optimal stabilizing control during vapor deposition is to select such source parameters as power q max , speed of his movement v(ξ), and shape parameter (heat flux width) H, for which the temperature of the silica tube is kept in the vicinity of a given state. We formulate the optimal control problem considering the linearization of the initial model and the further transition to the research of temperature deviations. Consider the problem of distributed control and a distributed monitoring system (11), (12) [19]. Let the function ∆q from the right side of the equation be the control function (further noted as ∆u), and observation is a function of the state of the system ∆Θ at every point in space and time. As a control space, we consider the space U = L 2 ((0, τ) × (0, L)), and as a space of solutions consider space Ω = (0, τ) × (0, L), where τ is the process control time. Notice that the value τ generally does not match the value T, defined above.
We define a target functional of an integral type that explicitly depends on the control function and has the property of coercivity, which will allow us to use the technique of obtaining an optimization system in the future and determine the optimal control function in the following explicit form: Functional F is a mixed integral criterion with an integrand that gives the average phase deviations and the total energy costs. Here is a positive numerical parameter (control price), which is estimated in advance or selected from the solutions of test problems.
Since control ∆u(t, z) is included in the task (11), (12) linearly, then we can talk about linear operator Λ acting on U, with values in the state space Ω, where Λ(∆u(t, z)) = ∆Θ(t, z)). Note that the Gato differential (the first variation of the mapping) on the increment δ∆u matches the value of the operator itself, due to linearity Λ, since We introduce the following notation: ∆ . Θ(t, z) = Λ(δ∆u); the value ∆ . Θ(t, z) will be called the conjugate state to the state of the system ∆Θ(t, z). We vary the initial differential problem (11), (12) and obtain a differential problem with respect to the conjugate state function of the following form: Multiply the right and left sides of the problem Equation (14) by an arbitrary function p(t, z) ∈ Ω and integrate over a given region to obtain the following: We obtain the following result after carrying out equivalent transformations: We require that the arbitrary p(t, z) function must satisfy the differential problem with the inverse time: Considering (16), the relation (15) transforms into the following: According to one of the consequences of the optimality criterion for distributed systems with compromise control [4], we have the following: where δF is the first functional variation (13), ∆u is minimizing element (optimal control function), and δ∆u is a variation of element ∆u.
Note that because of the notation introduced earlier, the integral of the first term in (18) exactly coincides with the integral on the left side of relation (17). With this in mind, Equation (18) can be rewritten as follows: Since, by definition, variation δ∆u = 0, the optimal control function can be represented by an explicit dependence as follows: Thus, considering the forms of differential problems (11), (12), (16), and the relation for the optimal control function (19), we can write the optimality system in strong form, i.e., in the form of a system of boundary value problems for the original function and the adjoint function in the following form:

General Results
The optimality system (20) is a differential problem that is inconsistent in time. The state of the system is particularly described by a boundary-value problem. It is correct in time forward, and the conjugate state is correct in backward time, which leads to some implementation difficulties: It is not possible to count in a certain direction in time. In this regard, solution (20) was carried out by the iterative method. In particular, knowing the values of the conjugate function p(t, z) at the i-th step, and hence the optimal control at the i-th step (i = 0, 1, 2, . . . ), we can solve the boundary-value problem that describes state of the system (in this case, we find the i-th approximation of the function ∆Θ(t, z)). The obtained solution allows us to solve the boundary value problem for finding (i + 1)th conjugate state approximations, where as a function ∆Θ(t, z), the i-th approximation is taken. The (i + 1)th optimal control approximation through the values of the conjugate function in the (i + 1)th approximation. A zero function approximation p(t, z) was chosen arbitrarily.
Each boundary value problem was solved by the finite element method implemented in the Comsol Multiphysics computing system. All physical functions and constants, which are the thermophysical parameters of the problem, were taken from the Material Library Comsol Multiphysics. To solve boundary value problems describing the state of the system and the conjugate state, a uniform grid in the spatial variable z was used. A study was conducted regarding the effect of the grid pitch on the resulting solutions. Moreover, when grinding the steps of the grid with respect to the spatial variable, starting from a certain value, the establishment of solutions was observed-the values of the desired functions changed very little (Figure 1). The dependence graph of the target functional on the iteration number, shown in Figure 2 is also interesting. It can be seen that starting from the seventh iteration, the functional decreases slightly and is set near the value F = 32,750 (for the value σ = 8 × 10 −8 ). It is assumed that the temperature distribution along the heating surface of the silica cylinder is controlled by a scanning pyrometer, which allows you to simultaneously record the temperature distribution over a length of 100 mm. Temperature measurements take place at regular intervals of time τ. The values of τ are selected considering the features of the MCVD process and physical capabilities to carry out the adjustment of the task parameters. Within the borders of this formulation, the system control time τ values were chosen equal to 1, 3, and 5 s. A corrective control was calculated according to Equation (19) from the solution of the system. Therefore, the heat source power q max is found, which is necessary to stabilize the process within the estimated time (0, τ). Note that the stabilizing control found will only act for the specified time period, after which the process repeats. A new temperature measurement takes place, the optimal control problem is solved again, the heat source power is corrected again, etc.
The scanning pyrometer takes readings in the heating zone only. As a demonstration example, we consider the difference between the real (determined using a pyrometer) and program temperature at a time corresponding to the coordinate z = L/2. We take the values of these differences as the initial conditions for the heat equation of the optimality system (20) (Figure 3).  Figure 4 shows the surface of the temperature distribution ∆Θ(t, z) and control ∆u opt (t, z) (heat flux) obtained as a result of solving the optimality system (20). The control price is the parameter, as noted above, that can vary. The values are selected from considerations of the method practical implementation. In our case, the choice of this value depends on the measurement error with a pyrometer. The used pyrometers have an instrumental error of 0.5%, so the control price parameter σ was selected on the assumption that the temperature measurement error is no more than 8 • C. Based on this, the problem was calculated for various values of the parameter σ ( Figure 5), and the green color shows the temperature profile, corresponding to the value of the error stated above.
Numerical experiments were conducted covering the influence of control time τ on the distribution of the optimal temperature (temperature ∆Θ(τ, z), corresponding to optimal control ∆u opt (t, z)). From the graphs presented in Figure 6, it is obvious that with increasing control time, the temperature values ∆Θ(t, z) are decreasing. The initial state of the system was chosen unchanged within the calculations.

Two-Dimensional (Axisymmetric) Mathematical Model of MCVD Process
We consider the MCVD optimal control problem in a two-dimensional (axisymmetric) formulation. The thermal state of the system, as before, will be described by Equation (2) from the beginning and boundary conditions (3) to (7). Then, we proceed to the twodimensional model linearization. The linearization technique described above will look like the following: where Θ(t, r, z) = Θ * (t, r, z) + ∆Θ(t, r, z) is with little variation ∆Θ(t, r, z) regarding distribution Θ * (t, r, z), which meets external influences q(t, z) caused by the operation of a heat source (1) on the boundary r = R.
If the function G is differentiable (at least twice) in the totality of arguments, then the linearization procedure leads to the following linear differential equation of the form concerning ∆Θ(t, r, z):

∂∆Θ(t,r,z) ∂z
We obtain the values of the corresponding derivatives: Based on the foregoing, the linearized equation takes the following form: Define functional operators Γ i , i = 1, . . . , 5, depending on the functions Θ(t, r, z), q(t, z), Θ 0 (r, z) and first derivatives of the function Θ(t, r, z) in spatial variables by the following relations: Similarly, we linearize the initial and boundary conditions: Then, the boundary and initial conditions of the linearized differential problem take the following form: So, the linear differential problem (22)-(27) is formulated, which is described by the two-dimensional axisymmetric heat equation.

Solution of the Boundary Optimal Control Problem
As noted above, the decisive role in the process of vapor deposition is played by the temperature distribution in the thermophoresis zone. In this regard, the most natural statement of the MCVD process optimal control problem is to obtain the temperature closest to the desired (specified) in the silica tube by controlling the burner power and/or speed.
In the scope of this statement for the optimal control function, we use the notation similar to the previous one-dimensional statement. The ∆u(t, z) optimal control function is renamed. In this case, it implies the power of the heat flux from the ∆q(t, z) heat source obtained by the system from the burner in contact with the outer surface of the silica tube (boundary r = R). The optimal control function ∆u(t, z) must be found from the condition of minimizing the functional of the integral form, under the condition that the positive parameter σ (control prices) is limited as follows: So, a functional of the form (28) minimizes temperature deviations from the programmed (specified) state on the outer cylindrical surface of the silica tube. At the same time, the control function is small in norm, and this means that the corrective actions on the system are small. The control efficiency in problem (22)-(28) consists in timely stabilization by adjusting the temperature regime supplied by the burner to the outer cylindrical surface of the silica tube. The control function here is the power of the heat flux from the burner, therefore, at the boundary r = R, the heat flux is set, and all types of heat exchange with the environment are considered. The system is also monitored on this part of the boundary; however, the parameter that is being observed is the temperature distribution on the surface of the cylinder. The stabilizing control problem formulated in this way is a problem with boundary control and boundary observation; the objective functional (28) for this problem, as before, clearly depends on the control function.
Based on the known variational principles and using the previously described algorithm for obtaining the necessary solvability conditions, as well as skipping a series of standard actions similar to those described earlier, we obtain an optimization system for the two-dimensional formulation of (22) to (28).
The optimality system is obtained in the strong form, i.e., in the form of a system of boundary value problems for the ∆Θ(t, r, z) original function and p(t, r, z) conjugate state function and has the following form:

∂∆Θ(t,r,z) ∂z
, ∆Θ(0, r, z) = ∆Θ 0 (r, z), λ ∂∆Θ(t,r 0 ,z) ∂r r λp(t, r, z) − ∂ ∂r (β 2 p(t, r, z))+ + ∂ 2 (λp(t,r,z)) The optimal control function is determined by the following formula: We assume that we know the desired temperature distribution in a cylindrical pipe in advance. This distribution can be obtained, for example, from experimental data or any other considerations. To obtain the temperature field closest to the desired in a silica tube, the following control algorithm is proposed. We measure the temperature field using a scanning pyrometer at equal intervals of time τ. As an initial condition for the problem of stabilizing optimal control (29), we consider the difference between the desired and actual temperatures. From the solution of the problem, we find the optimal value of the increment in the heat source power q max . This control will be valid only for the time interval (0, τ), after which a new temperature measurement takes place, etc. The objective of optimal control is to find a level of increase (decrease) in the power of the heat source at which the difference between the current temperature distribution and the program would be minimal. As a simulation object, we take a pipe made of synthetic silica with an external radius of R = 14 mm, an internal radius of r 0 = 12 mm, and a length of L = 1000 mm. Using a scanning pyrometer, we will measure the temperature. Figure 8 shows the profile of desired and actual (obtained by the measurement) temperature.
The difference in these temperatures will be the initial condition for the problem of stabilizing optimal control (Figure 9).
We decide the optimality system based on the obtained data. As a result of the solution, we obtain the optimal value of the heat source power increment q max. We will perform calculations using the COMSOL Multiphysics software product, as well as the MatLab and Maple software packages. To calculate the optimality system, we will use the simple iteration method with a constant control time τ and control cost σ. As a result of solving the optimality system (29), we obtained the temperature and power distribution of the burner over the time interval (0, τ).
Furthermore, from the obtained power function, it is necessary to obtain the optimal value of the power increment of the heat source q max . As noted earlier, the thermal flux of the MCVD flare is described by the Gauss function (1). In this regard, the burner physically cannot produce the heat flux described in Figure 10b, and the function ∆u must be approximated by a function that describes the shape of the heat source. As a result of approximation, we obtain the following law of change q max from the time (Figure 11).    Furthermore, the obtained values of the power source increment of the heat source must be converted to the corresponding values of the hydrogen flow rate at the MCVD installation. The resulting function will operate over time τ, after which a new temperature measurement will be made, etc. We also note that one of the key parameters affecting the control process is the control price σ. By changing this parameter, it is possible to accelerate the process of minimizing the module ∆Θ ( Figure 12).

Adjustment and Experimental Determination of the Parameters of a Moving Exposure Source
To determine the correspondence between the power of the heat source and the hydrogen flow rate, we carry out a series of experiments on heating silica pipes with a burner flame and measure using a scanning pyrometer at fixed intervals of the temperature distribution on the outer surface of the pipe. According to the received information, it is possible to select the desired values q max and H so that when they are substituted into Equations (22)-(27), the solution of the latter would give the temperature distribution on the outer surface of the pipe, which coincides with the experimentally obtained values. We assume that we know the law of motion of the source, i.e., the change in speed v(t) is known in advance (Equation (1)).
As noted earlier, to identify a moving heat source, it is necessary to determine the power function q max , form parameter H and source moving speed. We also know that when the flow rate of the burner hydrogen changes, the shape and power of the source changes. Figure 13 presents a comparison of simulation results and experimental data. As the flow rate of the hydrogen-oxygen mixture increases, the values of the parameters will also increase q max and H. The experiments showed that, in the researched range, the parameters q max and H are linearly dependent on hydrogen consumption (Figures 14  and 15, as well as Equations (31) and (32)).  So, it is easy to obtain the dependence of q max on dispersion (width) H. Returning to the optimal control problem, we obtain the distribution of the increment in the hydrogen flow rate from the process time ( Figure 16). We obtained a linear two-dimensional heat equation for the MCVD process. The problem of optimal stabilizing control with boundary control and boundary observation is formulated. An optimality system is obtained and an algorithm for its solution is proposed. The results of numerical studies of the temperature distribution ∆Θ and optimal control (burner power) ∆u opt are presented. The distribution of the increment of hydrogen flow from the control time is obtained.

Optimal Control Problem: A Simplified Mathematical Model of Fiber Drawing
The scheme of optical fiber drawing and its main geometric characteristics are presented in Figure 17. Due to the isothermal formulation of the problem, the solution area is a zone that is filled with a viscous Newtonian liquid. Let us descry a one-dimensional coordinate system (here, the spatial axis Oz is oriented downward). The direction of Oz is similar to the drawing direction (along the axis of symmetry). The technological drawing process can be described by the following system of partial differential equations (one-dimensional approximation): where z is the longitudinal coordinate; t is time; R(t, z) and V(t, z) are the silica melt flow velocity and radius, respectively; ρ and µ are the constant density and viscosity of silica melt, respectively; L is solution area length; R 0 is silica preform radius; R s (z) and V s (z) are functions corresponding to the initial state of the system; and V 0 and V L (t) are preform feed and fiber drawing velocities, respectively. Let it be that t ∈ [0, τ] and z ∈ [0, L]. Let us define by Ω the area within the segment [0, L], by Ω t Cartesian product [0, τ] × [0, L], by Γ t piecewise-smooth boundary of Ω t , here with Equations from (34) are satisfied inside the Ω t area, initial and boundary conditions from (34) are satisfied on parts of Γ t . Let us indicate the necessary spaces: L 2 (Ω t ) is the space of functions are summed squared on Ω t ; W 1 2 (Ω t ) is the space of functions from L 2 (Ω t ), which have a generalized first derivative belonging to L 2 (Ω t ).
Based on the state of system (34), let us represent the control problem. It is required to maintain the system near a state in which the deviations of the geometric characteristics of the elongated thread (fiber radius) from some optimal characteristics are minimal, provided that the fiber winding velocity is within the specified range. In this formulation of the problem under consideration, the control function will be the time function V L (t) defined above. At the same time, we will observe the radius of the finished fiber that is the function R(t, L) is the observed parameter.
Let us call the system movement at predetermined known values of velocities and radii as programmed movement, and the control corresponding to the movement as programmed control.
According to [20,21], if some nonlinear system is linearized in the neighborhood of its stationary state, then the analysis of such a linearized system can replace the analysis of the original system. In this case, the stationary value of the fiber radius is conveniently considered as the optimal geometric characteristic.
Let us linearize the system (34) in the following assumptions: where R(z) and V(z) are functions that correspond to programmed stationary movement; the function u(z) is a programmed control respective to the stationary mode; R(t, z) and V(t, z) are deviations of the factual movement from the programmed one; and u(t, z) is the deviation of the real control from the programmed one. Substituting the above relationships (35)-(37) into the system (34) and neglecting the nonlinear terms, let us obtain a linear problem in the following form: Here, the coefficient β(z) has the following form: The initial-boundary-value problem (38) obtained in this way describes the states that can be interpreted as deviations (disturbances) from the stationary programmed states of the original system (34).
Let us consider the functional on L + 2 (Γ 2 ) = {x(t) ∈ L 2 (0, τ) : x(t) ≥ 0} as follows: where α is a constant that is called the control price. Thus, the problem (38), (40) (with boundary control and boundary observation) consists in searching for the minimum value of the functional while the control parameter u takes values from the space of functions L + 2 (Γ 2 ) for α > 0. The optimality system obtaining as a result of linearization is the necessary conditions for the existence of the problem (38), (39). For this case, there is a possibility to obtain the optimality system in the boundary value problem form for the parent equations and equations for the conjugate states.
The questions of optimization problem solution existence for parabolic systems with objective functional in special forms was discussed in [4]. For the existence of at least one function u 0 (t) (optimal element) on which the functional (40) reaches its exact lower bound, it is necessary that the functional possess the properties of convexity, lower semicontinuity, and coercivity as follows R(t, z) = R(z)(1 + R(t, z)).
In accordance with the optimality criterion the value of the Gato differential [4] at the optimal element, u 0 equals to zero as follows: where . R = R u u= u 0 , · is the weak differentiation operator, and δ u 0 = w − u 0 is the control function variation. It should be noted that the differential problem in formulation (38) is linear with respect to the control function u(t). Thus, the state function of the system R(t, z) can be represented as a result of the action of some linear operator (note it as Λ) on the control function u(t) as follows: Using the linearity property of the operator Λ, let us obtain the following: Let us formulate the initial-boundary-value problem (38) [22] as follows: Let us multiply the continuity Equation (45) on a function q(t, z) ∈ L 2 (Ω t ) and the movement equation on a function p(t, z) ∈ L 2 (Ω t ), integrate both equations over the domain Ω t , and find their sum. The result of the addition will be the following expression: Let us use the Green's formula for terms containing derivatives of functions . R and . V as follows: Let us require that the functions q(t, z) and p(t, z) satisfy the following differential problem: Then, considering (45) and (46), the last integral relation takes the following form: Substituting it in (42): It follows that Finally, the optimality system takes the following form: The process of solving the system (47) was conditionally divided into some stages: 1.
Searching for a stationary solution of the system (34) (i.e., definition of functions R(z) and V(z)); 2.
Determination of a function β(z) depending on stationary states; 3.
Searching for a solution of the optimality system (47) and finding the optimal control function u 0 (t) (the numerical solution of the optimality system (47) was implemented using the Comsol Multiphysics modeling package); 4.
Result analysis.
Let us consider the process in detail. Under stationary conditions, the solution of the system of equations can be found in an explicit analytical form as follows: To solve the optimality system, the following values of the input parameters were used: τ = 5 s, ρ = 2200 kg/m 3 , µ = 10,000 Pa·s, L = 0.3 m. The function of the fiber radius deviation from its stationary state is shown in Figure 18 (line a). The deviation is specified in the form of an upward-convex function where the maximum is 10%. Special attention is paid to the choice of the control price value α. This parameter is necessary to determine in advance or based on the solutions of test problems. In this case, its value is equal to 0.5.
The optimal control function V(t, L) is shown in Figure 19. Let us solve the problem given above under the other initial condition. Let the function R s (z) be now the downward-convex function (line b in Figure 18), and the maximum of the fiber radius deviation from its stationary state is also 10%. Then, the optimal control function has the following form (line b in Figure 19).
To analyze the results, let us consider the obtained optimal control functions (lines a and b in Figure 19). As we can see, the optimal control functions are mirror symmetric with respect to each other, while the functions describing the initial state of the fiber radius are also mirror symmetric. It should be noted that the adjustment of the fiber winding speed is from −3 to +3 percent, which corresponds to the real production capabilities.
To estimate the influence of control on the resulted fiber quality, let us consider two drawing modes: Mode 1 is a non-control mode, and Mode 2 is a control mode according to the optimal state u 0 = 1 2α R − 3µ

Optimal Control Problem for the Optical Fiber Drawing in a Formulation That Takes into Account Surface Tension and Gravity Forces
Let us consider a system of partial differential equations that describes the process of drawing silica optical fibers in the following form: For this case, the system of equations describing the perturbation of the steady-state drawing mode (isothermal case) has the following form: where coefficients β 1 (z), β 2 (z), α 1 (z), and α 2 (z) have the following forms: , Here, ρ and µ are the density and viscosity of the melt, respectively, σ is the coefficient of surface tension, and g is the gravitational acceleration. The coefficients β 1 (z), β 2 (z), α 1 (z), and α 2 (z) are calculated through the functions R(z) and V(z), which are the solutions of the system of equations for the stationary problem (34).
Let us formulate the problem of boundary control and boundary observation of a parabolic system that is supplemented by corresponding initial and boundary conditions for the disturbed movement as follows: where the functions R(t, z), V(t, z), u(t) ∈ L 2 (Ω), The optimal control function u(t) must be found as a function that delivers a minimum to the integral functional under the condition that parameter α > 0 is small as follows: Furthermore, using the method described in the previous problem, let us obtain the optimality system in the following form: where q(t, z), p(t, z) ∈ L 2 (Ω) are given auxiliary functions and the optimal control function ∂p ∂z | z=L . The numerical implementation of the solution of the optimality system (52) was implemented in a similar way as described above:   Let us solve the same problem given above. Here, the initial deviation of the fiber radius R s (z) has the following form (line b in Figure 22), and the maximum of the deviation from its stationary state is also 10%. Then, having solved the optimality system (52) for the same values of the input parameters, we obtain the function V(t, L) in the following form (line b in Figure 23).
Let us analyze the results of solving the optimality system (52) with two different initial conditions for the equation of state. Note that the adjustment of the fiber winding speed is from −15 to +15 percent, which corresponds to the real production capabilities. To estimate the influence of control on the resulted fiber quality, let us consider two drawing modes: Mode 1 is a non-control mode, and mode 2 is a control mode according to the optimal state u 0 (t) = 1 ∂p ∂z | z=L . It is clear from Figure 24 that, in mode 2 (dotted line), the function of the resulted fiber radius stabilizes noticeably faster than in the mode without control (solid line). In addition, the deviations of the fiber radius are smaller in mode 2.

Conclusions
Overall, in this paper, two problems of optimal control of the processes of alloying of silica optical fiber preforms by the MCVD method and the silica optical fiber drawing are formulated, justified, and solved. The first problem was solved in two settings: onedimensional and axisymmetric, and the control function is the heat flux on the outer surface of the tube. In the second problem, two mathematical models are also considered: a simplified model and a model that takes into account the forces of surface tension and gravity. Both fiber drawing problems are solved in a one-dimensional formulation, where the control function is the winding speed of the resulted fiber. For both problems, linearized mathematical models were used, for which optimality systems were obtained in a strong form, i.e., in the form of partial differential boundary value problems. Each of the fiber drawing problems is solved for two different types of initial conditions that specify the deviation of the radius relative to its stationary (programmed) state. Numerical solutions of the optimality system, including the optimal control functions, are presented.