Optimal Design of a Delayed Feedback

In this paper robustness analysis and simultaneous stabilization of a semiconductor laser (SL) with a delayed optoelectronic feedback (OEF) is studied. Using the theory of stability radius maximal allowed norm bounded perturbations of the SL in contrast to stabilizable intervals of OEF parameters are obtained. Using proposed algorithms we can analyse the effect of interdependencies of uncertainties of the OEF and SL parameters and design an optimal OEF according to specific robustness constraints.


I. INTRODUCTION
Analysis of plants with time delays has attracted great interests in the last three decades [1].Time delays are present in many applications and areas [2], ranging from biology, process control, logistics, telecommunications, to semiconductor lasers (SLs) [3], which are indispensable in many sensor applications [3], [4] or optical communications [5], [6].
Especially when a SL is used in a specific application, its operating features, which depend on its stability, robustness and dynamics characteristics are of huge importance.Many characteristics can be affected by the use of OEF which can be applied to any laser diode [3].OEF is a connection from the optical laser output, which is detected by a photodetector to the injection current at the input.The detected photocurrent is amplified and fed back to the injected current as a positive or negative amount (Fig. 1).Typical result of such a feedback is shown in regular or irregular (chaotic) pulsations in the output.By appropriate setting of the feedback parameters stable and fast pulses with specific pulse width can be easily obtained [3].In this way high coherent light source can be achieved.Many applications use such a light source as for example laser range finders, markers, designators, night vision goggles, rifle sights, LIDAR, 3D optical data storage [7].In recent years the use of such systems has been attractive in secure chaos communications [6].The use of OEF introduces a delay to the system as a consequence of time responses of a detector and electronic circuits.The design and implementation of OEF is therefore certainly not trivial as different physical disturbances acting on closed loop system might affect the system to such an extent that stable output pulsation might become chaotic.
In this paper we are interested in the optimal design of the OEF in sense of robust stabilization and chaos suppression in the SL output [3], [8].We have followed the steps from [9] as far as computation of robustness of specific parameters is concerned and applied an improved method in the sense of efficient computation.We obtained maximal bounds of SL parameters uncertainties in contrast to simultaneously stabilizable intervals of OEF parameters.We propose numerical algorithms for derivation of such intervals inside of uncertainties sets.According to proposed algorithms and applied specific criteria optimal OEF design parameters can be determined, which guarantee certain amount of bounded perturbations of the SL and its OEF parameters.All the proposed algorithms were implemented using MATLAB ® and evaluated on a SL model [3].Presented results can be adapted and used for any other system that contains time delays.
The mathematical model of the SL with OEF is described in Section II.Stability and robustness analysis is presented in Sections III and IV.Simultaneous stabilization and optimal OEF design is addressed in Section V. Design examples and conclusion are given in Sections VI and VII.

II. SEMICONDUCTOR LASER WITH OPTOELECTRONIC FEEDBACK
The dynamics of a SL with OEF (Fig. 1) can be described with two equations of the photon number and the carrier density.The mathematical model is represented with the following rate equations [3]: where  is the photon number,   is the gain coefficient,  is the carrier density,  ℎ is the carrier density at threshold,  is the injected current density,  is the thickness of the active layer,  is the elementary positive charge constant,  is the feedback strength,   is the constant offset in the feedback loop,   is the steady state value for the photon number,  0 is the carrier density at transparency,   is the lifetime of the carrier, and  is the feedback time delay.See [10] for the details of those parameters.Local dynamics of SL with OEF might be investigated using linear stability analysis.According to [3], [8] the linearization procedure applied to (1) and ( 2) results in the following characteristic equation in variable  with one which can be represented as a delay differential equation (DDE) [11]   where  =  −0.5 ,  is the ratio of the photon damping rate in the cavity to the rate of population relaxation,  =  0 /(√1 + ), where  0 = √ − 1, where  = (  )/(    ) is the pumping rate,   is steady state carrier density and  = /(  ) is the feedback coefficient.Values of the parameters [3] are listed in Table I.From the stability analysis of (3) boundary of stable and unstable oscillations of the laser can be determined (Fig. 2).
As was reported in [8] this curves correspond to the stable and unstable boundaries of the rate equations (1) and ( 2).As will be shown in our further analysis local stability boundaries might rapidly change by perturbations of certain parameters.Beyond these bounds the laser output is subjected to various irregular oscillations, which extend from regular pulsing over quasi-periodic pulsing to chaotic pulsing, where chaotic pulses have both chaotic peak intensities and chaotic pulse intervals.To assure stable output laser oscillations despite certain perturbations of parameters analysis of robustness must be taken into account.

III. STABILITY ANALYSIS
Stability of a time delay system can be determined by its eigenvalues, which should be located in the open left half plane [2], [12].There are several methods that provide numerical algorithms for eigenvalues computation of a time delay system [13]- [18].Stability boundaries concerning feedback gain in contrast to time delay were determined using an improved version of [14] according to the Algorithm I.
The [14] is based on the computation of all eigenvalues in the selected half plane via spectral discretization and on the estimation procedure of the desired discretization points.In addition, a correction method of derived eigenvalues is used.But there is a drawback, the region of the computation must be set manually.The larger is the selected region the more expensive and time consuming is the computation.On the other hand, the computation is unsuccessful if the region doesn't contain any eigenvalues.For determination of the stability region of (3) only the rightmost eigenvalue is needed.We have improved the algorithm [14]  3. Iteratively select  = ( 1 +  2 )/2 using bisection algorithm until   <  <   is met.4. Calculate eigenvalues using spectral discretization method [14].
Using the Algorithm I stability margins for negative   and positive   gain values of (3) were determined for the following set of time delays (Fig. 2 and Fig. 3) IV. ROBUSTNESS ANALYSIS Perturbations of certain parameters might rapidly change stability bounds and influence regular SL output oscillation to such an extent that it might become chaotic.In what follows, maximal allowed perturbations (ℎ) [9] were computed for stable region of OEF parameters (  ,  , ) using proposed implementation.

A. Stability Radius
Maximal allowed perturbations of a time delay system, which shift its eigenvalues to imaginary axes, can be determined using complex stability radius, which is defined by [9]     where   is a structured singular value and can be computed using the mussv routine in MATLAB ® , matrices  and  define certain structure of perturbations,  −   − ∑    −   =1 for  =  is a characteristic equation of a time delay system with system matrices   , … ,   .Applying additive perturbations to (4), where where  ∈ [0, 1], results in following matrices:

B. Efficient Computation of Stability Radii for the Semiconductor Laser
Computation of numerous stability radii on a 2-D grid (Fig. 3) of stabile region is a demanding task due to the complex nature of structured singular value problem, which in addition needs to be evaluated over an interval of  (6).The computation was performed over the set (5) and following stable sets of feedback gains To reduce the complexity of computation we propose Algorithm II. ) −1 ‖ 2 on an interval of  1 = [0,10/  ] with step size  1 = 0.01/  .2. Extract two thirds of the largest  1 with corresponding  1 as  2 and  2 .3. Compute stability radius  1 (6) over the  2 with step size of  1 .4. When the algorithm stops, rectangle with the largest  ∆ is selected and saved, as well as   and   .Step size was constructed as a function of time delay to assure appropriate precision also for higher values of time delays.With higher delays eigenvalues condense towards the origin of the complex plane [2], which results also in condensed picks in magnitude.For that reason smaller step size is needed.
Figure 3 shows results of the stability radii for the stable region (Fig. 2).The highest levels of stability radii are displaced from the center of the stable region.With increased values of the stability radii stability margins start to narrow until individual connected contours are derived.

V. SIMULTANEOUS STABILIZATION
Another important aspect of robustness concerns the time delay.As a result of different external influences the time delay in the OEF might change and cause instability.Figure 3 clearly denotes the reduction of the areas surrounded by individual contours when the amounts of perturbations increase.The higher are perturbations of parameters the more limited become stability bounds and the admissible intervals of gains and delays decrease.

A. Simultaneous Stabilization
Computation of maximal stabilizable intervals ( ∆ ,  ∆ ) of feedback parameters (, ) in contrast to maximal allowed parameters perturbations (ℎ(, )) was carried out for all contours of stability radii inside of regions of  1 and  2 (Fig. 2).Individual regions are defined as: , The goal was to obtain the largest possible rectangular shape inside of each contour (Fig. 3 and Fig. 4).By choosing the centred points (  ,   ) of such rectangles optimal feedback parameters (  ,   ) can be obtained which allow bounded parameters perturbations (ℎ) and simultaneous stabilization on delay ( ∆ ) and gain intervals ( ∆ ) of certain desired size.
To obtain   , ,   , and corresponding  ∆ , for each contour   and for  ∆ , we propose Algorithm III.
Algorithm III.Computation of    ,    and  ∆  intervals for  ∆  = {0.005|= 1,2, … }: 1.Each contour   = {0.005|= 1,2, … } of ℎ , was extracted, up sampled to delay and gain step size of 0.0001 and divided into upper and lower set of points.2. For each point  the corresponding ,  and  points were found (Fig. 4). was found using a bisection algorithm.Figure 5 and Fig. 6 show results of the Algorithm III evaluated over contours in region  1 .Figure 5 represents inverse dependence of  ∆ and   to  ∆ .
Figure 6 represents   ,   for which  ∆ and  ∆ at specific   can be achieved.For higher values of   smaller stabilizing intervals can be achieved, which clearly shows that it is impossible to obtain maximal robustness of the SL and OEF.

B. Interpolation of Families of Curves
Resulting curves of (  ,   ) and ( ∆ ,   ) can be represented by several groups of polynomials, which can be derived using curve fitting tool cftools in MATLAB ® .
Family of curves of (  ,   ), belonging to region  1 (Fig. 6), bounded by the area   ∈  1 ∩  2 ∩   > 0 ∩   ≥ 0 , where  1 and  2 are derived using fourth order polynomial fitting, is described as a group of polynomials with polynomial coefficients defined as . ∆ ( ∆ , ) belonging to  1 , bounded by  ∆ ∈  3 ∩  > 0.07 ∩  ∆ > 0.01, where  3 is derived using linear line fitting, is described using polynomials In order to achieve accurate fitting we can divide the whole family of curves into several groups of polynomials (12) or (14).Polynomials ( 12) and ( 14) were derived using proposed Algorithm IV.
Algorithm IV.Interpolation of a family of curves: 1. Family of curves is divided into several groups according to  ∆ ,  or  ∆ ,   .2. For each group interpolate all curves and save coefficients.
3. Interpolate evolution of individual coefficients of polynomials to derive (13) and describe the whole group with ( 12) or ( 14) respectively.

C. Correction Algorithm
Obtained polynomials   ( ∆ ,   ) were used to correct corresponding polynomials  ∆ ( ∆ , ) with the help of the Algorithm I and proposed Algorithm V. Correction of each polynomial was executed for 15 distinct points.
parameters in contrast to uncertainties of the SL parameters, which was evaluated using several algorithms.Improved numerical algorithms for computation of stability radius and stability boundaries in terms of computation of only few rightmost eigenvalues were provided.As was shown, amounts of simultaneously stabilizable OEF intervals inversely depend on the values of SL parameters perturbations, which imply that it is impossible to obtain maximal robustness of the SL and the OEF parameters.In this manner an optimal OEF might be designed considering specific constraints and criteria as discussed in provided example.
Proposed algorithms can be adapted and used for any other time delay system.
in a manner to derive only few rightmost eigenvalues by means of automatic selection of computation region.In this way an efficient algorithm for determination of stability boundaries of (3) has been constructed in MATLAB ® .Algorithm I. Computation of few rightmost eigenvalues: 1. Initialize computation region {() ≥ | ∈ ℂ},  = −/3 and execute estimation procedure of discretization points  of [14].If   = 10 <  <   = 100 is met go to step 4. 2. If  >   iteratively decrease region  =  + 5 and save  1 =  until   <  <   or if  <   save  2 = .If  <   iteratively increase region  =  − 5 and save  1 =  until   <  <   or if  >   save  2 = .If   <  <   is met go to step 4.

3 .
The formation of rectangles proceeds until points occur inside of the rectangle.Only first time rectangles are calculated for all upper points.Next iterations use the information of the position of previously selected maximal rectangle in Step 4. 4. When the algorithm stops, rectangle with the largest  ∆ is selected and saved, as well as   and   .

Fig. 4 .
Fig. 4. Algorithm III: Out of blue rectangles the red rectangle is chosen with maximal interval of delays.