A Direct Method for Determining a P-Solution of Linear Parametric Systems

For given functions (1b), the iterative method of [11] is applicable if the radius of p is not larger than a threshold ra (the applicability radius of the method considered [10]) within which the iterations are convergent. A shortcoming of the method of [11] for determining the p-solution x(p) is the fact that the number of iterations, needed to compute it, can be prohibitively large for relatively wide parameter vectors p whose r(p) is close to ra. In an attempt to improve its computational efficiency, a direct method for determining the p-solution is suggested in the present paper. Its functional characteristics are first compared with those of a direct method for computing an outer interval (nonparametric) solution of LIP systems (1). The new direct method is also compared with the known iterative method of [11].

A new type of solution x(p) to the LIP system (1) (called parameterized or p-solution) has been recently introduced in [11]. It is defined as a corresponding linear interval form where L is a real n×m matrix while l is an n-dimensional interval vector. An iterative method for determining x(p) was suggested in [11].
As shown in [11], the p-solution x(p) is an enclosure of S 1 , i.e.
That is why, the use of x(p) seems rather promising in solving various global optimization problems where the LIP system (1) is involved as an equality constraint.
For given functions (1b), the iterative method of [11] is applicable if the radius of p is not larger than a threshold r a (the applicability radius of the method considered [10]) within which the iterations are convergent. A shortcoming of the method of [11] for determining the p-solution x(p) is the fact that the number of iterations, needed to compute it, can be prohibitively large for relatively wide parameter vectors p whose r(p) is close to r a . In an attempt to improve its computational efficiency, a direct method for determining the p-solution is suggested in the present paper. Its functional characteristics are first compared with those of a direct method for computing an outer interval (nonparametric) solution of LIP systems (1). The new direct method is also compared with the known iterative method of [11].

Preliminaries
We will use boldface to denote interval quantities, underscores to denote lower bounds and overscores to denote upper bounds. Subscripts will be used to denote components of vectors or matrices. In general, vectors (scalars) will be denoted by lower case letters, while matrices will be denoted by upper case. Let A and B be n×m real matrices. Relations A=B, A≤B etc. are meant component-wise. The same convention is accepted for interval matrices. An interval n×m matrix A can be defined either in its lower-upper-end form The inverse H of an interval matrix A, denoted symbolically H=A -1 , is the hull of the set G={A -1 : A∈A}. The special case of an interval n×m matrix where I is the identity matrix will be used in the sequel. Let where σ(∆) is the spectral radius of ∆. Let The following result is known. The inequality (5) is a necessary and sufficient condition for B to be regular [12].
The LIP system (1) can be written equivalently in the form( where Å, A (µ) , µ=1,…, m are n×n real matrices, A 0 is a n×m real matrix, p is a symmetric interval vector ( 0) =  p and ˆ(1,..,1) = T p is an m-dimensional unit vector. The collection of all real matrices A(p) when p varies within p will be denoted A(p) and will be referred to as interval parametric (IP) matrix. A similar notation b(p) will be also used for the collection of b(p), p∈p. Thus, (7a) can be written symbolically as The narrowest interval matrix containing A(p) will be called (interval) hull of A(p) and will be denoted A(p). Thus Obviously The Direct Method

Basic version
The direct method for determining x(p) is based on Lemma 1 and comprises the following steps.
First, we assume (temporarily) that the IP matrix A(p) is regular so 1 − =  R A exists. Using R, (7a) is transformed equivalently as in [6]: From (9)ˆ( C is the µ th column of C 0 ). It is seen from (7), (10) that the original problem of finding the p-solution x(p) has been reduced to finding the p-solution y(p) of (10). Indeed, if the united set solution of (10) is denoted S 2 =∑(B(p), b(p), p), then from (7), (9) and (10) Now an interval (nonparametric) matrix B=B(p) of the type (4) is introduced using (10a): Next, (10) is enclosed in p by the following "mixed type" system It is assumed that so B is regular. Let S 3 denote the united solution set of (13). Obviously,

B(p)∈B for any p∈p so
Indeed, the parameter dependence between the left-hand and right-hand side of (10) is ignored in (13). However, (13) is much easier to handle than (10).
and is replaced with Due to (14), the matrix H can be computed using (6) where H  is a diagonal matrix whose non-zero elements are all positive In view of (18), equation (17) is rewritten in the form and it is seen that L is real n×m matrix while s=[-s, s] is a symmetric interval vector.
Finally, it will be shown that the p-solution sought is given by (ii) the p-solution x(p) of the given LIP system (1) exists and is determined by (20), (21).

Proof:
The assertion (i) follows from the regularity of B, the inclusion B(p)∈B for all p∈p and the relation B(p)=RA(p). The assertion (ii) is proved as follows. From (11) and (15) The implementation of the present direct method based of Theorem 1 will be referred to as algorithm A1. As is easily seen, the bulk of the computation is related to computing the m matrices B (µ) so the number of arithmetic operations (multiplications) N 1 is approximately

Improved version
A better version of the direct method is possible which is based on the concept of the applicability radius r a of an interval method [10]. Following the general case approach [10], we introduce the family of parameter vectors of variable width where p 0 is given (start) vector and ρ is a variable scalar. Obviously, the direct method is applicable as long as the IP matrix A(p(ρ)) remains strongly regular. Let B(ρ) denote the interval matrix in (12) as a function of , i. e. B(ρ)=I+∆(ρ)[-1,1]. Hence, On account of (14) As is well known, σ(∆(ρ)) is defined as the eigenvalue λ m of maximum magnitude from the eigenvalue problem It is readily seen from (12) that ∆(ρ)=ρ∆. Now for ρ=ρ 0 =1, But B(ρ) becomes singular for so r a =ρ 1 and from (25) Using the concept of applicability radius, we have the following result.

Theorem 2: Let 
A in (7b) be nonsingular. Let the applicability radius r a of the direct method considered have been computed by (26). Then, for any ρ< r a : (ii) the p-solution of system (1) associated with p∈p(ρ) exists and can be determined by (20), (21).
In practice, it is sometimes necessary to apply algorithm A1 repeatedly for various input vectors p(ρ) according to (22). In such cases, a better version of A1 is possible which is based on Theorem 2 and the relations (ρ)=ρ∆, B 0 (ρ)=ρB 0 .
Algorithm A2: It comprises two stages.  x(p(ρ)) for each p(ρ) will be found using

Comparison with Other Methods
We compare the present method (referred to as method M3) with the iterative method of [11] (method M2) and the direct method of [6] (method M1) according to the following three criteria: a) enclosure efficiency: tightness of the approximation of the solution set S 1 obtained by the respective p-solutions , b) computational efficiency, c) applicability radius. The algorithms of methods were programmed in MATLAB environment. The programs were run on a 1.7 GHz PC computer.

Comparison with the direct method [6].
We first compare M2 with M1 whic is also a direct method. It should, however, be stressed that method M1 yields an outer solution of (1) in the standard (nonparametric) form of an interval vector x. We now prove the following result which is a corollary of Theorem 1.

Corollary 1:
The hull x(p) of the p-solution of x(p) of (1) obtained by method M3 is equal to the outer solution x of (1) obtained by method M1, i.e.

x(p)=x. (27)
Proof: It suffices to prove that v(p)=v. According to [6] (using the present paper notation) is the so-called companion matrix of B. Since p and b are symmetric, According to the present direct method, formula (20) and the fact that It is seen from (28) and (29) that r 3 =r 1 which completes the proof.
In proving Corollary 1, we have also shown that the assumption, adopted in [6], for B to be an H-matrix is superfluous.
The advantage of method M3 over M1 reveals itself when solving problems other than finding an outer solution x of (1). (A vast class of such more general problems, referred to as the generalized interval hull solutions, has recently been defined in [10].) To illustrate the above assertion, we consider the following parametric linear programming (PLP) problem (the simplest possible representative of this class) where the constraint is the LIP system A(p)x=a(p) with [11] The parameter interval vector p is given by its centre and radius p 0 =(0.5 0.5 0.5), r 0 =(0.5 0.5 0.5). (30c) Using p 0 and r 0 , we transform (30b) into the equivalent form (7b) (as shown in [11]) to have p µ ∈[-1,1]: where 0.5 1.5 -0.5 1.5 -3 0.5 , (1 0.5 1) 1.5 3 1 The range of (30a) is the interval  (A(p), a(p), c(p), p)

={f=c T (p)x: A(p)x=a(p), p∈p}
(denoted for shortness as f * ). For simplicity, we have chosen (in the general case, c=c(p) and c(p) can be nonlinear functions). What we seek is to find an outer bound f on f * .
In the case of M1, f is given as where x i are the components of the interval outer solution of (1). The bound f obtained by M3 and denoted f 3 is found as the range of Unlike (32), f 3 is determined by (33a) using the components x i (p) of the p-solution x(p) of (1).
To show quantitatively that (33b) is narrower than (32) we employ the merit figure where r(f 1 ) and r(f 3 ) are the radii of f 1 and f 3 , respectively. On account of Corollary 1 hence, from (32) Thus, To find r(f 3 ), we write x i (p) as From (36) and (37) It is seen that method M3 has better enclosure efficiency than method M1.
Using (22), we now compute η 31 for various ρ within the applicability radius r a of method M3. The radius r a is determined as follows. Solving the partial eigenvalue problem (24) for ρ=1, we have found the corresponding 0 1.3425 λ = m . By (26) The corresponding values for η 31 are given in the second row of Table 1.
It is seen that η 31 decreases as ρ grows. This is explained by the fact that the relative weight of the first term ∑ j j L in (37) decreases with respect to the second term s 0 in function of ρ.

Comparison with the iterative method [11]
According to [11] x(p) is obtained iteratively by computing linear interval forms l (k) =c (k) +L (k) p+s (k) at each k th iteration. These forms enclose corresponding solution sets ( 1) 2 k S + . It is shown that ( 1) 2 k S + tends to S 2 as k grows to infinity (if the iteration process is convergent). So encloses S 1 of (1). This approach takes into account the interdependency between all the parameter components p µ of p. The present direct method only accounts for the parametric dependencies in the righthand side b(p) in (10). Therefore, the method of [11] is expected to be better according to criterion a).
On the other hand, the computation volume of the direct method is much smaller than that of the iterative method since the amount of computation needed in the former method is roughly the same as that required on each iteration of the latter method. Thus, the iterative method is bound to be more expressive than the direct method and this discrepancy will become more pronounced as the radius of p approaches the applicability radius r a . Let x (2) (p) and x (3) (p) denote the p-solutions of (31) obtained by methods M2 and M3, respectively. The enclosure efficiency of the two p-solutions will be compared solving the PLP problem (30a), (31). In this case, we compare the outer solutions f 3 (given by (33)) with f 2 determined in a similar way using x (3) To quantitatively assess the superiority of M2 over M3, we use the merit figure To show the dependence of η 23 on the parameter width ρ, we need the applicability radius r a (M2)  (the iterative process becomes divergent for ρ=0.72). It is seen from (39) and (43) that the present direct method M3 has a larger applicability radius than the iterative method M2.
We now show the dependence of η 23 on ρ from ρ=0.1 up to ρ=0.7.
The corresponding values of η 23 are given in the second row of Table 2. As expected, the iterative method M2 provides tighter outer bound f (2) (p) as compared to f (3) .
The two methods are also assessed as regards their computational efficiency using the index '  Table 2 show that the better enclosure efficiency of method M2 is obtained at the cost of much large computer times, especially when using algorithm A2, for ρ close to r a (M2). Also r a (M2)<r a (M3). Therefore, in some cases, it may be preferable to use M3, A2 rather than M2.

Conclusion
A direct method (method M3) for determining the p-solution of the LIP system (1) has been suggested in the present paper (Theorem 1). Using the concept of applicablyity radius r a , a better version of the method has been proposed (Theorem 2). It is proved (Corollary 1) that the hull x(p) of the p-solution x(p) of (1) obtained by method M3 is equal to the outer solution x of (1) obtained by the method of [6] (method M1). Method M3 is, however, superior to method M1 in solving problems other than finding an outer solution x of (1).As an illustration, the parametric linear programming (PLP) problem (30a), (30b), (31c) is considered. The numerical results obtained show that M3 provides tighter PLP solutions than M1 (Table 1). Method M3 is also compared with the iterative method of [11] (method M2) ( Table  2) using the PLP problem. It has been shown that M3 has a larger applicability radius r a (M3)> r a (M2), is much faster than M2 (especially when the second algorithm of M3 is employed), providing p-solutions x (3) (p) that are wider but comparable in width with the solutions x (2) (p) of M2. For these reasons, it may be preferable, in some cases, to use M3 rather than M2.
In view of the theoretical considerations and the numerical evidence in the paper, it is expected that the use of the p-solution x (3) of LIP systems (1) can lead to the development of new more efficient methods for solving various global optimization problems [10,11].
Appendix: A procedure for computing (I-ρ∆) -1 for an arbitrary ρ<r a is suggested here. It is based on the relation between M=I-∆ and the corresponding inverse P=M -1 , on the one hand, and M ′ =I-ρ∆ and its inverse (I-ρ∆) -1 , on the other. We first transform where ( 1) : and ( 1) : are the ith column and row of P (i-1) , respectively. Finally, on account of (A1) to (A4) P ρ =[ρ(M+αI)] -1 =(1/ρ)P (n) . (A5) As can be easily seen, updating P (i-1) to P (i) costs n 2 multiplications so the total modification of P (0) to P (n) requires N=n 3 multiplications. This is a better result as compared with inverting M ′ by a standard method.