A parallel algorithm for semi-infinite programming

https://doi.org/10.1016/S0167-9473(03)00152-XGet rights and content

Abstract

The implementation results of an algorithm designed to solve semi-infinite programming problems are reported. Due to its computational intensity, parallelisation is used in two stages of the algorithm—for evaluating the global optimum and for checking the feasibility of constraints. The algorithms are parallelised using MPI—the message passing interface. The algorithms are then applied to engineering and macroeconomic modelling problems, including designing robust optimally informative experiments for dynamic model identification and verification, and robust macroeconomic policies for inflation targeting.

Introduction

A parallel implementation of the solution of the following semi-infinite problem is considered:minx∈Xf(x)s.t.G(x,y)⩽0,∀y∈YX={x∈Rn|gi(x)⩽0,i=1,…,k}Y={y∈Rm|qi(y)⩽0,i=1,…,l}G:Rn+mRnG,where X and Y are nonempty compact subsets of Rn and Rm, nG is the number of coupled constraints Gj(x,y)(j=1,…,nG) and f and G are twice continuously differentiable on X and X×Y, respectively. The term semi-infinite programming arises from the fact that G(x,y)⩽0,∀y∈Y represents an infinite set of constraints on x.

Problem (1) was first considered by John (1948) who gave necessary and sufficient conditions for its solution. Since then, it has been extensively studied in the literature. Other efforts include Blankenship and Falk (1976), Coope and Watson (1985), Polak and Mayne (1976), Mayne and Polak (1982), Gustafson (1981) and, more recently, Lawrence and Tits (1998). The reader is referred to Hettich and Kortanek (1993) for further references.

Problem (1) can be used to formulate certain minimax and saddle point problems arising in optimal engineering design under parametric uncertainty. Applications to computer-aided design and parametric optimal decisions are presented in Kwak and Haug (1976) and Mayne et al. (1979). In such problems, f(x) may represent the cost function where x is a vector of control/design variables and y is a vector of uncertain parameters. The formulation ensures that the resulting design will fulfil the constraints G(x,y)⩽0 for all yY, usually related to the inherent safety of the system. Worst-case design is essentially a minimax formulation when the best decision is taken in view of the worst-case scenario. This can be formulated as a semi-infinite program. It has the advantage that the decision is robust to uncertainty. Robustness is assured as system performance is optimal for the worst-case scenario and will be improved if any other scenario is realised (Rustem and Howe, 2001; Zakovic and Rustem, to appear).

In this paper, (1) is solved using the algorithm based on the work of Blankenship and Falk (1976) and Darlington et al. (1999), where a two-phase procedure to guarantee feasibility for the general nonlinear constraints was employed. The algorithm described in Darlington et al. (1999) is specialised towards constraint satisfaction in view of the worst-case realization of the uncertainties. This is further discussed in (Zakovic and Rustem, to appear) where the problem is solved using a two stage procedure that searches for global maximum violation of the constraints, together with a version of the algorithm that searches for any violation of constraints.

The main difficulty of the algorithm is the computational effort required. This is mostly in the computation of the global optimisation and the checking of constraint violation. In the present study, a parallel scheme is provided to address these difficulties and the advantages are illustrated using numerical examples.

The algorithm begins at k=0 with a finite subset YkY, solving the approximate problem:vk=minx∈Fkf(x),Fk={x∈X|G(x,y)⩽0,∀y∈Yk}.F and v are defined asF={x∈X|G(x,y)⩽0,∀y∈Y},v=minx∈Ff(x).The algorithm solves the kth auxiliary global optimisation problem—compute ymaxY such thatG(xk,ymax)⩾G(xk,y),∀y∈Y,where xk=argminx∈Fkf(x). Termination is declared when the inequalityG(xk,ymax)⩽0is satisfied. Otherwise, the algorithm setsYk+1=Yk+{ymax},k=k+1,and solves (2) again.

Note thatF⊂⋯⊂Fk+1⊂Fk⊂⋯⊂F0,impliesv⩾⋯⩾vk+1⩾vk⩾⋯⩾v0.To justify the stopping rule, observe that the solution of the kth auxiliary problem satisfiesG(xk,ymax)⩽0,thenG(xk,y)⩽0,∀y∈Y.Thereforevk⩽v⩽f(xk)=vk,so that xk is a solution of (2), with Fk=F.

Algorithms for global optimisation problems are very complex in nature as they usually involve a large number of iterations, functions evaluations and local searches. Dixon and Szego 1975, Dixon and Szego 1978 distinguish two basic approaches—deterministic and stochastic. Most of the methods to solve the global optimisation problem can be placed in these two categories.

Deterministic methods find the global minimum by locating all the minima of the objective function. The global minimum is then the observed local minimum with the lowest function value. Methods like grid search (space covering), trajectory methods or interval arithmetic methods (i.e Hansen, 1992; Wolfe, 1996) are classified as deterministic methods, because, for a given initial point, they always find the same solution. Such methods usually involve additional assumptions on f. Interval arithmetic methods can in principle solve the global optimisation problem within some desired error bound (Hansen, 1992; Moore, 1979). Methods like random search (Deep and Evans, 1994) and clustering (Rinnooy Kan and Timmer, 1987) are classified as stochastic, as they are based on random sampling in the region of interest and given an initial point may find a different solution if the run is repeated.

Recent years, especially as the required computer power has become available, have seen rapid development developments in global optimisation algorithms. Most notably among these are, to name only a few are those of Floudas et al. Adjiman 1996, Adjiman 1998, Smith and Pantelides 1997a, Smith and Pantelides 1997b, or a paper on parallel search global optimisation algorithms by Pardalos (1989). A good survey on recent development and trends in the area can be found in Pardalos et al. (2000).

For global optimisation in (3), a stochastic algorithm, developed by Rinnooy et al. (1987) is used. It is a well known algorithm and has been used, among other things, for solving complex real life problems arising in laser diffractometry (Zakovic, 1997; Zakovic et al., 1998)

The algorithm starts independent local optimisations from numerous different points, and keeps a list (Xmax) of maxima encountered. As a stopping rule, it uses a Bayesian estimate of the total number of maxima, based on the numbers of local optimisations performed, and local solutions found. Also xmax(xst)=argmaxx{f(x)} in step 4 is defined as a local maximum of f(x), obtained when the local optimisation was started from xst. Algorithm 1 summarises Rinnooy Kan and Timmer's method.

In step 10 there is a possibility that nloc=nmax+2. In such the case ntot is not computed, but another local optimisation is performed. If new maximum is not encountered then nlocnmax+2, and ntot is well defined. However if there is a new maximum, meaning that still nloc=nmax+2, the whole procedure is repeated until there has been performed a local optimisation which did not encounter new maximum.

The condition f(xmax(xst))⩾0 in step 11 is not part of the original Rinnooy Kan and Timmer algorithm. The price of computing the global optimum can be, and usually is, very high. As the search for points which violate constraints is performed, it may well be computationally cheaper to stop the global optimisation search at any point where violation is encountered. This approach can bring benefits regarding the computational time, especially when functions that are being globally optimised are very complex and nonconvex Zakovic and Rustem (to appear). There may be the setback, however, that by including local violations of constraints, the nonlinear programming problem in (2) is solved more times than necessary. A switch is built into the proposed algorithm, which allows the user to decide, depending on the complexity of the functions f(x) and G(x,y), which version of the algorithm to use.

The constrained nonlinear programming problems that appear in (2), and in Algorithm 1, are solved using standard nonlinear programming packages; in this case the NAG subroutine e04ucf.2 The subroutine computes the minimum of a general smooth function that may include simple bounds as well as linear and nonlinear constraints.

Pseudo code for solving (1) is summarised in Algorithm 2.

The feasibility check in step 11 depends on the nature of the constraints, whether they are hard or soft, and on the application being considered. Hard constraints indicate situations where no constraint violation can be tolerated by the physical system. This requires feasibility for all realizations of y, which may further increase the price of computation. If strict feasibility cannot be obtained or guaranteed, a penalty approach may be used, allowing optimal solutions with the smallest possible violation of the constraints G(x,y) within Y. This approach is further discussed in Darlington et al. (1999).

As can be seen, the algorithms used to solve the semi-infinite programming problem are computationally intensive. In order to accelerate their solution, parallelisation techniques are used; in particular, the proposed algorithms are parallelised using the Message Passing Interface (MPI). MPI is a specification for a library of routines to be called from C or Fortran programs, and was first developed in 1993 when a group of computer vendors, software writers, computational scientists, and others collaborated on setting a standard portable message passing library.

Parallel computers have recently become an everyday tool of computational scientists, and have widespread applications in optimisation and mathematical programming (Pardalos et al., 1992; Pardalos, 1999; Migdalas et al., 1997). However, there are still problems that prevent the widespread use of parallelism, both in terms of hardware and software. For instance, for hardware, it is still difficult to build intercommunication networks that can keep up with the speed of most advanced single processors. For software, compilers that automatically parallelise sequential algorithms are very limited in their applicability.

Parallelisation brings up some issues that are, obviously, not present in a serial context. Most important among those issues are task allocation, where the breakdown of the total workload into smaller tasks is performed so that such tasks are allocated to different processors. Another issue is communication, of intermediate results of different processors. Also among these important issues is synchronisation of the different computations processes.

MPI's design and development has been guided not only by the above issues but by the need for portability and to impose no semantic restrictions on efficiency. It has also been designed to be a convenient and complete definition of the message passing model (Gropp et al., 1994). The message passing model consists of processes that have only local memory, but are able to communicate with other processes by sending and receiving messages. In the message passing model, data transfer from the local memory of one process to the local memory of another process requires operations to be performed by both processes.

Although it cannot be said that the message passing model is superior to any other parallel methods, it has become widely used for a number of reasons, such as universality, expressivity, ease of debugging, performance, to name but few (Gropp et al., 1994). Efficiency is one reason too, for instance, some Crays have hardware barriers. Although there are about 125 functions in the MPI standard, there are only 6 functions that are indispensable, upon which a large number of efficient programs can be written. All the other functions add flexibility, robustness or convenience.

Fig. 1 illustrates task farming which is a communication pattern that can be implemented using message passing. The master process, enumerated Process 0, sends data to the slave processes, enumerated 1 to n−1. After the data is processed it is sent back to the master process.

There are two sub-procedures in the full SIP algorithm that can be parallelised, namely the global optimisation in step 10, of Algorithm 2 and checking for violation of the constraints in step 11.

It is worth mentioning that stochastic approaches are becoming very important tools in a wide range of fields, starting from physics, chemistry or applied mathematics to social sciences. Bearing this in mind it is clear that an important aspect of such approaches is the source of pseudo random numbers. A scalable package for parallel pseudo random number generation (SPRNG) is being developed at Florida State University.3 In this implementation however, a generator from NAG library (see footnote 2) is used, as the number of starting points generated is relatively small.

In the case of global optimisation, parallelisation is implemented as follows—the master process randomly chooses Nst starting points xst and subsequently distributes them evenly (≈Nst/n−1 points) to each of the remaining n−1 processors. Local optimisations are performed by each processor, each starting from one of the points in the starting point set xst. The results are then sent back to the master process, where it is checked whether the newly obtained locally optimal points are disregarded or put in the recorded list of local minima.

A similar procedure is followed for the constraint violation problem. There are nG constraints present, where each can be sent to a different process, which subsequently returns the information about the violation. Then, again, the master process decides whether to add a new constraint to the feasible region.

Section snippets

Applications

To illustrate the application of the algorithms discussed in this paper, two examples are presented, one from macroeconomics and the other from engineering, detailing and discussing their solutions.

Conclusions

A parallelised semi-infinite programming algorithm has been presented in this paper. Two different sub-procedures—global optimisation and the checking of constraint violation have been parallelised. The parallelisation has been effected within the MPI (Message Passing Interface) framework (Gropp et al., 1994). The hardware used is Fujitsu AP3000 parallel server—a distributed memory parallel server based on powerful 64-bit UltraSPARC technology. The nodes of the AP3000 are connected by Fujitsu's

References (38)

  • J.W. Blankenship et al.

    Infinitely constrained optimization problems

    JOTA

    (1976)
  • I.D. Coope et al.

    A projected Lagrangian Algorithm for semi-infinite programming

    Math. Programming

    (1985)
  • Deep, K., Evans, D.J., 1994. A parallel random search global optimisation method, Computer Studies 882, Loughborough...
  • Dixon, L.C.W., Szego, G.P. (Eds.), 1975. Towards Global Optimisation, North Holland,...
  • Dixon, L.C.W., Szego, G.P. (Eds.), 1978. Towards Global Optimisation, Vol. 2, North Holland,...
  • W. Gropp et al.

    Using MPI

    (1994)
  • S.A. Gustafson

    A three-phase algorithm for semi-infinite programming

  • E. Hansen

    Global Optimisation Using Interval Analysis

    (1992)
  • L. Hansen et al.

    Robust permanent income and pricing

    Rev. Econom. Stud.

    (1999)
  • Cited by (6)

    • Automatic scenario generation for efficient solution of robust optimal control problems

      2024, International Journal of Robust and Nonlinear Control
    • Relaxation-based bounds for semi-infinite programs

      2008, SIAM Journal on Optimization
    1

    Present address: Orbis Investment Advisory Ltd, 5 Mansfield St, London W1G 9NG.

    View full text