Mathematical programming package HYBRID

Marek Makowski IIASA, Laxenburg, Austria"

Janusz S. Sosnowski

Systems Research Institute, Polish Academy of Sciences, Warsaw.

Abstract

HYBRID is a mathematical programming package which includes all the func- tions necessary for the solution of multicriteria LP problems and single-criteria linear-quadratic problems. HYBRID is specially useful for dynamic problems since the applied algorithm exploits the structure of a dynamic problem and the user has the advantage of handling a problem as a dynamic one which results in an easy way of formulation of criteria and of interpretation of results. HYBRID is ori- ented towards an interactive mode of operation in which a sequence of problems is to be solved under varying conditions (e.g., different objective functions, reference points, values of constraints or bounds). Criteria for muitiobjective problems may be easily defined and updated with the help of the package. Besides that HYBRID offers many options useful for diagnostic and verification of a problem being solved.

HYBRID is available in two versions: one for VAX 6210 (running under Ultrix-32) and one for a PC compatible with PC IBM! AT!XT.

1 Introduction

The purpose of this report is to provide a sufficient understanding of mathematical, methodological and theoretical foundations of the HYBRID package. Section 1 con- tains executive summary, short program description and general remarks on solution techniques and package implementation. Section 2 contains mathematical formulation of various types of problems that can be solved by HYBRID. Section 3 presents method- ological problems related to solution techniques. Section 4 presents foundations of the chosen solution technique and documents the computational algorithm. Section 5 con- tains short discussion of testing examples. Last two sections contain conclusions and references.

'on leave from the Systems Research Institute of the Polish Academy of Sciences, Warsaw.

This paper does not include information necessary for using the package. A reader who is interested in usage of the package should consult a User Guide to HYBRID (Makowski and Sosnowski, 1988b). In the User Guide the following topics are discussed:

• the way of choosing various options provided by the package.

• guidelines for formulation and modification of a problem which is to be solved or at least processed by HYBRID.

• the way in which HYBRID provides diagnostics and results.

• a short tutorial example.

• the specification of the MPS standard for input data and an example of the MPS format input file.

1.1

Executive summary

HYBRID is a mathematical programming package which includes all the functions necessary for the solution of linear programming problems. The current version of HYBRID, called HYBRID 3.1, may be used for solving both static and dynamic LP problems (in fact also for problems with a more general structure then the classical for- mulation of dynamic linear problems). HYBRID 3.1 may be used for both single- and multi-criteria LP problems as well as for single-criteria linear-quadratic problems. Since HYBRID is designed for real-life problems, it offers many options useful for diagnostic and verification of a problem being solved.

HYBRID is a member of the DIDAS family decision analysis and support sys- tems since it is designed to support usage of multicriteria reference point optimization.

HYBRID can be used by an analyst or by a team composed of a decision maker and an analyst or-on last stage of application-by a decision maker alone. In any case, we will speak further on about a user of HYBRID package.

HYBRID can serve as a tool which helps to choose a decision in a complex situation in which many options may and should be examined. Such problems occur in many situations, such as problems of economic planning and analysis, many technological or engineering design problems, problems of environmental control. To illustrate possi- ble range of applications, let us list problems for which the proposed approach either has been or may be applied: planning of agriculture production policy in a decentral- ized economy (both for governmental agency and for production units, Makowski and Sosnowski, 1985a), flood control in a watershed (Kreglewski et aI., 1985), planning for- mation and utilization of water resources in an agricultural region, scheduling irrigation, planning and design of purification plant system for water or air pollution.

To avoid a possible misleading conclusion that the usage of HYBRID may replace a real decision maker, we should stress that HYBRID is designed to help the decision maker to concentrate on his actual decision tasks while HYBRID takes care on cumber- some computations and provides information that serves for analysis of consequences of different options or alternatives. A user is expected to define various alternatives or

scenarios, changing his preferences and priorities When learning about consequences of possible decisions.

HYBRID could be used for that purpose as a "stand alone" package, however- after a possible modification of a problem in an interactive way---one can also output the MPS-format file from HYBRID to be used in other packages. The later approach can be used also for a transformation of a multicriteria problem to an equivalent single- criteria LP. HYBRID includes also some diagnostic functions that are not performed by many other linear programming packages, e.g., by MINOS (it is interesting to note that the authors of MINOS actually advise the user to debug and verify the problem with another package before using MINOS).

HYBRID can be used for solving any linear programming problem but it is spe- cially useful for dynamic problems; this covers a wide area of applications of operation researches. Many optimization problems in economic planning over time, production scheduling, inventory, transportation, control dynamic systems can be formulated as linear dynamic problems (Propoi, 1976). Such problems are also called multistage or staircase linear programming problems (Fourer, 1982, Ho and Hanne, 1974). A dynamic problem can be formulated as an equivalent large static LP and any commercial LP code may be used for solving it, if the problem corresponds to single objective optimization.

For multicriteria problems, a preprocessor may be used for transformation of that prob- lem to an equivalent LP one. One of the first versions of the system DIDAS was a package composed of a preprocessor and a postprocessor for handling transformation of multicriteria problem and processing results respectively (Lewandowski and Grauer, 1982). Those pre- and postprocesors were linked with an LP package. HYBRID 3.1 has generally a similar structure. The main difference is that-instead of an LP package-- another non-simplex algorithm is applied, which exploits the dynamics of a problem and that HYBRID is integrated package with user-friendly interface. Similarly as some other systems of DIDAS family, HYBRID has the advantage of handling a problem as a dynamic one which results in an easy way of formulation of criteria and of interpre- tation of results, since one may refer to one variable trajectory contrary to a "static"

formulation of dynamic problems which involves separate variables for each time period.

HYBRID has been designed more for real-world problems that require scenario anal- ysis than for academic (e.g., randomly generated) problems. Thus HYBRID is oriented towards an interactive mode of operation in which a sequence of problems is to be solved under varying conditions (e.g., different objective functions, reference points, values of constraints or bounds). Criteria for multiobjective problems may be easily defined and updated with the help of the package.

The binary files with HYBRID 3.1 are available from IIASA in two versions: one for VAX 6210 (running under Ultrix-32 ver. 3.0) and one for a PC compatible with IBM/AT/XT.

1.2 Short program description

1.2.1 Preparation of a problem formulation

A problem to be solved should be defined as a mathematical programming model.

Formulation of a mathematical programming model is a complex task and this paper

is not devoted to discuss this question in detail. Therefore this section is aimed at providing only a short summary of a recommended approach.

Firstly, a set of variables that sufficiently describe the problem-for the sake of the desired analysis-should be selected. It is desired-however not necessary-to define the model in such a way as to possibly exploit the problem structure (further on referred to as a dynamic problem). Secondly, a set of constraints which defines a set of admissible (i.e. acceptable or recognized as feasible by a decision maker) solutions should be defined.

Finally a set of criteria which could serve for a selection of a solution should be defined.

The formal definition of criteria can be performed in HYBRID in an easy way.

However, it should be stressed that any definition of a complex model usually requires cooperation of a specialist-who knows the nature and background of the problem to be solved-with a system analyst who can advise on a suitable way of formal definition. It should be clearly pointed out that a proper definition can substantially improve the use of any computational technique. For small problems used for illustration of the method, it is fairly easy to define a model. But for real life problems, this stage requires a close cooperation between a decision maker and a team of analysts as well as a substantial amount of time and resources.

For real life problems, the following steps are recommended:

1. Mathematical formulation of the problem being solved should be defined.

2. A data base for the problem should be created. This may be done on PC with a help of a suitable commercial product (such as Framework, dBase, Paradox, Oracle or Symphony). Original data should be placed in this data base. A user need not worry about possible range of quantities (which usually has an impact on computational problems) because HYBRID provides automatic scaling of the model.

3. Verification of the data base and of the model formal definition should be per- formed.

4. The corresponding MPS standard file should be created. This may be done by a specialized model generator (easily written by a system analyst), or an universal generator such as GEMINI (developed at IIASA), or GAMMA (part of FMPS package on UNIVAC), or LPL (cf Hurlimann, 1988), or by any appropriate utility program of data base software. We strongly discourage the user from creating the MPS file with help of a standard text editor.

1.2.2 Model verification

This stage serves for the verification of model definition which is crucial for real appli- cation of any mathematical programming approach.

First stage consists of preprocessing the MPS file by HYBRID, which offers many options helpful for that task. HYBRID points to possible sources of inconsistency in model definition. Since this information is self-explaining, details are not discussed here. Itis also advisable to examine the model printout by rows and by columns, which

helps to verify model specification and may help in tracing possible errors in MPS file generation.

Second stage consist of solving optimization problems for selected criteria which helps in the analysis of consistency of solutions. For larger problems, the design and application of a problem oriented report writer is recommended. HYBRID generates a

"userJile" for that purpose which contains all information necessary for the analysis of a solution.

After an analysis of a solution, a user may change any of the following parameters:

values of coefficients, values of constraints and also any parameters discussed in next section. This may be done with help of the interactive procedure which instead of MPS file uses "communication region" that contains problem formulation processed by HYBRID. Therefore, a user needs no longer to care about original MPS file which has the backup function only.

1.2.3 Multiobjective problem analysis

For a given model, the user can define various multiobjective problems to be analyzed.

Problem analysis consist of consecutive stages:

• analysis of obtained solution

• modification of the problem

• solution of modified problem.

Analysis ofa solution consists of following steps (some of which are optional):

1. The user should examine of values of selected criteria. Since the solution obtained in HYBRID is Pareto optimal, the user should not expect improvement in any criteria without worsening some other criteria. But values of each criterion can be mutually compared. It is also possible to compute the best solutions for each criterion separately. A point (in criteria space) composed of best solutions is called the "utopia" point (since usually it is not attainable). HYBRID provides also a point composed of worst values for each criterion. This point is called

"nadir" point. Such information help to define a reference point (desired values of criteria) because it is reasonable to expect values of each criterion to lie between utopia and nadir point.

2. The user may also make at this stage modifications to the original problem without involving the MPS file.

3. For dynamic problems, HYBRID allows also for easy examination of trajectories (referred to by so called generic name of a variable).

Modification of the problem may be done in two ways:

1. At this stage, the user can modify the formulation of the original model. But main activity in this stage is expected after the model is well defined and verified and no longer requires changes in parameters that define the set of admissible (acceptable) solutions. Itshould be stressed, that each change of this set usually results in change of the set of Pareto-optimal solutions and both utopia and nadir points should be computed again.

2. If the values of all constraints and coefficients that define the admissible set of solutions are accepted, the user should start with computations of utopia point.

This can be easily done in an interactive way. After utopia and corresponding nadir points are obtained (which requires n solutions of the problem, where n is the number of criteria defined) the user can also interactively change any number of the following parameters that define the selection of an efficient solution to the multicriteria problem:

- Reference point (Le. desired values for each criterion) might be changed. This point may be attainable or non-attainable (cfsect. 2.4).

- Weights attached to each criterion can be modified.

Reference trajectories in dynamic case can be changed as reference points.

- Regularization parameters in selection function can be adjusted.

3. Additionally, the user can temporarily remove a criterion (or a number of criteria) from analysis. This option results in the computation of a Pareto optimal point in respect to remaining "active" criteria, but values of criteria that are not active are also available for review.

Solution

0/

a problem. The multiobjective analysis problem defined by a user (after possible modification) is transformed by HYBRID to an equivalent LP problem which is solved without interaction of a user (an experienced user may however have an access to the information that characterizes the optimization run).

1.2.4 Remarks relevant to dynamic problems

HYBRID allows for solving both static and dynamic LP models. Static models can be interpreted as models for which a specific structure is not recognized nor exploited.

But many real life problems have specific structure which-if exploited-can result not only in much faster execution of optimization runs but also remarkably help in problem definition and interpretation of results.

Numerous problems have dynamic nature and it is natural to take advantage of its proper definition. HYBRID offers many options for dynamic models, such as:

1. In many situations, the user may deal with generic names of variables. A generic name consists of 6 first characters of a name while 2 last characters corresponds to the period of time. Therefore, the user may for example refer to the entire trajectory (by generic name) or to value of a variable for a specific time period (by full name). Such approach corresponds to a widely used practice of generating trajectories for dynamic models.

2. The user may select any of 4 types of criteria that correspond to practical ap- plications. Those can be defined for each time period (together with additional

"global" conditions), but this requires rather large effort. Therefore, for dynamic problems, criteria are specified just by the type of criterion and the generic name of the corresponding variable. Types of criteria are discussed in detail later.

3. A model can be declared as a dynamic one by the definition of periods of time.

For a dynamic model, additional rules must be observed. These rules correspond to the way in which the MPS file has to be sorted and to the way in which names for rows and columns are selected. These rules follow a widely accepted standard of generation of dynamic models. The formulation of a dynamic model, which is accepted by HYBRID, is actually an extension of the classical formulation of a dynamic model (cfSection 2.2.). In our formulation, a model may contain also a group of constraints that do not follow the standard of state equations.

1.2.5 General description of the software package and data structure The package is constructed in modules to provide a reasonably high level of flexibility and efficiency. This is crucial for a rational use of computer resources and for planned extensions of the package and possible modification of the algorithm.

The package consists of five subpackages:

• Two preprocessors that serve to process data, enable a modification of the model, perform diagnostics and may supply information useful for the verification of a model. The first preprocessor is used for processing of initial formulation and di- agnostics of the model. It also transforms a multicriteria problem to a parametric single criteria optimization problem. The second preprocessor allows for analysis of a solution and for the interactive change of various parameters that may corre- spond to choice of some option, change of parameters in definition of multicriteria problem, change of matrix coefficients, right hand sides of constraints etc.

• Optimization package called solver of a relevant optimization problem (either static or dynamic).

• Postprocessor that provides results in the standard MPS format and generates the "user file" which contains all information needed for the analysis of a solution;

the later option makes it easier to link HYBRID to a specialized report-writer or a graphic package.

• Driver, which eases the usage of all subpackages. The PC version of driver provides a context sensitive help which helps an inexperienced user in efficient usage of the package.

All five subpackages use a binary file that contains all data defining the problem being solved. A second binary file contains a solution obtained by last run of the solver.

From the user point of view, HYBRID 3.1 is still one package that may be easily used for different purposes chosen via specification file.

The chosen method of allocating storage in the memory takes maximal advantage of the available computer memory and of the features of typical real-world problems. In general, the matrix of constraints is large and sparse, while the number of all non-zero coefficients that take different numerical values is much smaller than the number of all non-zero coefficients. A super-sparse-matrix technique is therefore applied to store the data that define the problem to be solved. This involves the construction of a table of coefficients which take different numerical values. The memory management is handled by a flexible way. HYBRID is coded partly in C and partly in an extension of Fortran (the latter part is processed by a preprocessor to generate a code which conforms to Fortran 77 standard). Such approach results in a faster (much faster for PC version running under DOS) execution and in a decrease of memory requirements.

Special commands of HYBRID support model verification and problem modification.

This is necessary to facilitate scenario analysis and to reduce the problems caused by inappropriate scaling (cfsect. 4.7).

The data format for the input of MPS file and the output of LP results follows standards adopted by most commercial mathematical programming systems (cf e.g.

Murtagh, 1981, Makowski and Sosnowski, 1988b).

1.2.6 Outline of the solution technique

HYBRID uses a non-simplex algorithm - a particular implementation of the augmented Lagrangian (or Lagrange multiplier) method - for solving linear programming prob- lems. General linear constraints are included within an augmented Lagrangian function.

The LP problem is solved by minimizing a sequence of quadratic functions subject to simple constraints (lower and upper bounds). This minimization is achieved by the use of a method which combines the conjugate gradient method and an active constraints strategy.

In recent years many methods oriented for solving dynamic linear problems (DLP) have been developed. Most of those methods consists of adaptation of the simplex method for problems with a special structure of constraints. In HYBRID, a different approach is applied. A DLP, which should be defined together with a state equation, is solved through the use of adjoint equations and by reduction of gradients to control subspaces (more exactly, to a subspace of independent variables). The method exploits the sparseness of the matrix structure. The simple constraints (lower and upper bounds for non-slack variables) for control variables are not violated during optimization and the resulting sequence of multipliers is feasible for the dual problem. The global con- straints (i.e. constraints other then those defined as simple constraints) may be violated, however, and therefore the algorithm can be started from any point that satisfies the simple constraints.

The solution technique can be also used to solve single-criteria quadratic problems with virtually no changes in the algorithm. However, a routine to input and handle the relevant data and a corresponding standard for data input have yet to be designed and implemented. So far only single criteria linear-quadratic problems in the form discussed in Section 2.5 may be solved. The solution method for multi-criteria quadratic problems requires modification of the algorithm. However the necessary modifications will be

based on HYBRID 3.1.

1.2.7 General description of options provided by the package

In order to provide general information about capabilities of HYBRID, the main options are listed below. HYBRID offers the following features:

• Input of data and the formulation of an LP problem follow the MPS standard.

Additional rules (that concern only sequencing of some rows and columns) should be observed in order to take advantage of the structure of a dynamic problem. An experienced user may speed up computations by setting certain options and/or parameters (cfthe HYBRID User Manual).

• The problem can be modified at any stage of its solution (i.e., by changing the matrix of coefficients, introducing or altering right-hand sides, ranges or bounds).

• The multicriteria problem is formulated and solved as a sequence of paramet- ric optimization problems modified in interactive way upon analysis of previous results.

• The solution technique can be chosen. First choice is done by definition of a static or a dynamic problem. Some specialized techniques may be used for badly con- ditioned problems that usually cause numerical problems. This includes one of two regularization techniques (see Section 4.5) and/or possibility of using precon- ditioned conjugate gradient method (cfSection 4.6). For a badly scaled problem, an implementation of scaling algorithm is available (as described by Makowski and Sosnowski (1981) and briefly discussed in Section 4.7).

• Comprehensive diagnostics is implemented, including the checking of parallel rows, the detection of columns and rows which are empty or contain only one entry, the splitting of columns, the recognition of inconsistencies in right-hand sides, ranges and bounds, and various other features that are useful in debugging the problem formulation. The package supports a display of a matrix by rows (printing the nonzero elements and names of the corresponding columns, right-hand sides and ranges), as well as a display of a matrix by columns (analogous to displaying by rows). A check of the feasibility of a problem prior to its optimization is optionally performed. More detailed information for an infeasible or unbounded problem is optionally provided by the package.

• All data that correspond to the formulation of the problem being solved are stored in a binary file. An other binary file contains all other information corresponding to a current run. The latter file is stored on disk in certain situations to allow continuation of computations from failed (or interrupted) runs or to run a mod- ified problem while using previously obtained information. Therefore the MPS input file is read and processed only by first preprocessor, which serves for ini- tial formulation of the problem. Such approach allows also for efficient storing of many solutions that may be later used for more detailed analysis, comparisons and modifications.

• Any solution is available in the standard MPS format and in a binary file which contains all data that might be useful for postoptimal analysis and reports.

1.3 Remarks on implementation

HYBRID 3.1 is an extended version of HYBRID 3.03 documented in (Makowski and Sosnowski, 1988a, 1988b). Therefore there are only small changes in the methodolog- ical guide in comparison to the methodology presented in (Makowski and Sosnowski, 1988a), because the solution techniques are basically the same. However, there are some important methodological innovations. The main differences are the following:

• The code has been modified as to allow for solution of single criteria linear- quadratic problems.

• The preconditioned conjugate gradient technique for minimizing augmented La- grangian has been implemented.

• The second regularization option which allows for finding the optimal solution with minimum distance from a given reference point has been made operational.

• The optimization algorithm has been improved by an automatic evaluation of some parameters, a different technical implementation of scaling, some changes in control flow, which results in its faster execution.

• The user interface (for PC version of the code) has been improved. A new approach to usage of the package and to data handling provides for easier use of the package.

• Diagnostics have been improved and several observed bugs have been removed.

• Part of the code has been rewritten in C language. This allows for more efficient memory management and usage. Change in the way of internal data handling resulted in remarkable improvement of execution speed.

2 Statement of optimization problems

2.1 Formulation of an LP problem

We will consider a linear programming problem (P) in the following standard form (see, e.g., Murtagh and Sanders, 1977):

minex b - r ~ Ax ~ b

l~x~u

(1) (2) (3) where x,e,l,u E R", b,rE Rffl and A is an m x n matrix.

The constraints are divided into two groups: general constraints (2) and simple constraints (3). In the input data file (MPS file) the vectors bis called RHS and the vector r -RANGES. The vector land u are called LOWER and UPPER BOUNDS,

respectively. Obviously, some of bounds and/or ranges may have an infinite value.

Therefore HYBRID may be used for solving any LP problem formulated in the way accepted by most of commercial packages.

2.2 Classical formulation of a dynamic LP problem (CDLP)

Before discussing a formulation of a dynamic problem that can be solved by HYBRID 3.1, let us first consider a classical formulation of a dynamic linear programming problem (CDLP) (cf Propoi, 1976) in the following form:

Find a control trajectory

U=(U ... ,UT) and a state trajectory

x

=

(Xl, ... , XT) satisfying the state equations with initial conditionXo

and constraints

(4)

FTxT

:s

dT which minimize the performance index

t=l, ... ,T

t

=

1, ... ,T (5)

(6) (7)

T

l]atxt

+

btUt)

t=1

where:

- t = 1, ... ,T denote periods of time

(8)

- state variables Xt, control variables Ue, both for each period, are elements of Eu- clidian spaces of appropriate dimensions;

- matrices At,Bt,Dt , Ft are assumed to be given,

- RHS vectors Ct and de, as well as range vectorrt and bounds for control variables et and It are given,

- initial condition Xo is given.

The above given formulation has been chosen for the purpose of simplification of presentation only. Actually, the following modifications are accepted:

1. Instead of inequality (5), equality constraints can be used;

2. Since no constraints of bounds type (6) are allowed for state variables x, such constraints may be specified in columns section of MPS file, thus formally are handled as inequality constraints of type (5);

3. Performance index (goal function) can either be specified as single objective or will be replaced by a dummy goal function that is defined by the transformation of a multicriteria problem to a parametric LP problem;

The structure of an CDLP problem (formulated above as in Propoi, 1976) may be illustrated by the following diagram (example for T

=

^{3, Uh}U2, U3, XO,XI,X2, X3 are vectors, slack variables are not shown):

Ul U2 U3 Xo Xl X2 X3 rhs var.

B l 0 0 Ao -1 0 0 Cl state eq.

0 B 2 0 0 Al -1 0 C2 state eq.

0 0 B 3 0 0 A 2 -1 C3 state eq.

D l 0 0 Fo 0 0 0 do constr.

0 D 2 0 0 Fl 0 0 d_l constr.

0 0 D3 0 0 F2 0 d2 constr.

0 0 0 0 0 0 F3 d3 final state

b_l b2 b3 0 al a2 a3

-

^goal

where1 is identity matrix and 0 is a matrix composed of zero elements.

2.3 Formulation of a dynamic problem (DLP)

The formulation of CDLP has been chosen for the purpose of simplification of presen- tation only. Actually HYBRID 3.1 is capable to solve problems of more general class, which will be referred to as Dynamic Linear Programming problems (DLP). Namely, the matrices B

=

diag(B,), D

=

diag(D,}, F

=

diag(F,) need no longer be block diagonal matrices. Also matrices below identity matrices need no longer have any specific struc- ture. Therefore the CDLP is a specific example of DLP. One of main generalizations- from a practical point of view-is that a problem with delays for control variables (which is not CDLP-class problem) may be solved by HYBRID. In fact, HYBRID accepts also problems with delays for both state and control variables, provided that state variables for periods "before" initial state do not enter state equations. A choice of criteria for CDLP-class problem is also limited in comparison with that for DLP (efsect. 4.3).

All variables are divided into two groups: decision variables^U and state variables Xt, the latter are specified for each period of time.

A single criteria DLP problem may be formulated as follows:

Find a trajectory Xt and decision variables ^U such that both:

state equations:

t - l

-Htxt

+ L

At-l"x,

+

BtU

=

et, ,=0

t

=

1, ...,T (9)

with given initial condition Xo

and constraints:

T

d - r ~

L

^Ftxt

+

^Du ~ d t=o

e~u~f

are satisfied and the following function is minimized:

(10) (11)

(12) Components of vector u are called decision variables for historical reasons. Actually a vector u may be composed of any variables, some of them may be specified for each time period and enter criteria defined for a dynamic case. But some components of vector u may not be specified for any time period (cfsect. 7.3.1). An example of such variable is ".. dummy.", a variable generated by HYBRID for a multicriteria problem.

A user may also specify variables independent of time. For the sake of keeping the formulation of the problem as simple as possible we have not introduced a separate name for such variables.

The following two symbols can be used in the specification file for definition of DLP:

NT - number of periods (stands for T in the above formulation)

NSTV - number of state variables in each period (the dimension of vectors Xt)

The user can define state inequalities instead of state equations (9). The slack variables for such inequalities are generated by HYBRID. Therefore, for the sake of the presentation simplicity, only the state equation will be considered further on.

The structure of an DLP problem may be illustrated by the following diagram:

(corresponding to an example analogous to the above example for CDLP)

u Xo Xl X2 X3 rhs var.

BI A_oo -HI 0 0 CI state eq.

B2

A

lO All -H₂ 0 C2 state eq.

B3 A20 A21 A₂₂ -H3 C3 state eq.

D Fa FI F₂ F₃ d constr.

b 0 al a2 a3 - goal

where H_t is diagonal matrix and 0 is a matrix composed of zero elements.

2.4 Multicriteria optimization

2.4.1 General remarks

The specification of a single-objective function, which adequately reflects preferences of a model user is perhaps the major unresolved difficulty in solving many practical

problems as a relevant optimization problem. This issue is even more difficult in the case of collective decision making. Multiobjective optimization approaches make this problem less difficult, particularly if they allow for an interactive redefinition of the problem.

The method adopted in HYBRID 3.1 is the reference point approach introduced by Wierzbicki (1980). Since the method has been described in a series of papers and reports and has been applied to DIDAS (ef Kallio et aI., 1980, Lewandowski and Grauer, 1982), we give only general outline of the approach applied. This approach may be summarized in form of following stages:

1. The user of the model (referred to further as the decision maker-DM) specifies a number of criteria (objectives). For static LP problem a criterion is a linear combination of variables. For DLP problems one may also use other types of criteria (cfsect. 2.4.2). The definition of criteria in HYBRID can be performed in an easy way described in the User Manual.

2. The DM specifies an aspiration level ii

=

{iil" ..^,iiNC},where iii are desired values for each criterion and NC is a number of criteria. Aspiration level is called also a reference point.

3. The problem is transformed into an auxiliary parametric LP (or DLP) problem.

Its solution gives a Pareto-optimal point. Ifspecified aspiration level ii is not at- tainable, then the Pareto-optimal point is the nearest (in the sense of a Chebyshev weighted norm) to the aspiration level. Ifthe aspiration level is attainable, then the Pareto-optimal point is uniformly better then ii. Properties of the Pareto- optimal point depend on the localization of the reference point (aspiration level) and on weights associated with criteria.

4. The DM explores various Pareto-optimal points by changing either the aspiration leveliior/and weights attached to criteria or/and other parameters related to the definition of the multicriteria problem.

5. The procedure described in points 3 and 4 is repeated until satisfactory solution is found.

To give more formal presentation, let us introduce following notation:

N C is the number of criteria qi is the i-th criterion

iii is the aspiration level for i-th criterion

Wi is a weight associated with i-th criterion (whereas the user specifies its absolute value which is internally changed to negative depending on the type of criteria--ef sect. 2.4.3).

em is a given non-negative parameter.

A Pareto-optimal solution can be found by the minimization of the achievement scalarizing function in the form

NC

. max (Wi(qi - iii))

+

~m

L

^{Wiqi -+ min}

.=I,... ,NC i=1

This form of achievement function is a slight modification of a form suggested by A. Lewandowski (1982) and by A. Wierzbicki (1978). Note that for ~m

=

0 only weakly Pareto-optimal points can be guaranteed as minimal points of this function. Therefore, the use of a very small ~m results (except of situations in which reference point has some specific properties) in properly Pareto-optimal solution with trade-off coefficients bounded approximately by _~mNC and _1/~mNC. If~m is very small, these properly efficient solutions might practically not differ from weakly efficient (Pareto optimal).

On the other hand, too big values of~mcould drastically change properties associated with the first part of the scalarizing function.

2.4.2 Types of criteria

A user may define any number of criteria. To facilitate the definition 6 types of criteria are available and a user is requested to declare chosen types of criteria before their actual definition. Two types of criteria are simple linear combination of variables and those criteria may be used for both static and dynamic problems. Four other types of criteria correspond to various possible performance indices often used for dynamic problems.

Since the latter criteria implicitly relate to the dynamic nature of the problem, they may be used only for variables that are defined for each time period. The only exception is the type DER of criteria, which may be defined by state variables only.

For the sake of simplicity, only the variables of the typeXi (which otherwise is used in this paper to distinguish a state variable in DLP) are used in the following formulae.

Note that Xi

=

^{XiI}' t

=

^{1, ...},T.

An k-th criterionqk is defined in one of following ways, for static and dynamic LP:

Type MIN

T n

qk

= L L

^{ail XiI -+ min}

1=1 i=1

where n is number of (state and control) variables, T is number of periods; T = 1 is assumed for static LP.

Type MAX

T n

qk

= L L

^{ailXil -+ max}

1=1 i=1

The following four criteria types are exclusively for dynamic LP:

Type SUP

qk = max (XiI - :tie) -+ min

1=1, ... ,T

where Xi is a selected state or control variable, Xi -its reference trajectory Type INF

q/<

=

min (Xit - Xie) ----. max l=l ...T

Type FOL

q/<= max (abs(xi/ - XiI)) ----. min

l=l, ... ,T

Type DER (which applies only to state variables)

q/< = max (abs(xil - xiI-d) ----. min

l=l ...T

2.4.3 Transformation of multicriteria problem to an auxiliary LP

The transformation is done by HYBRID 3.1, therefore its description here has only informative purpose. This description may be useful in case of using the MPS file (optionally created after modifications and transformation of a problem) as input for another LP package.

Following notation is used throughout this subsection:

v - name of the auxiliary variable v

Wi - optional weight coefficient for i-th criterion (default value equal to1.), eni - name of i-th criterion,

chI - string (2-characters) which identifies t-th period of time, iii - reference point (aspiration level) for i-th criterion,

qi - linear combination of variables that defines a criterion of the type MAX or MIN, , , - delimiters of a string,

T - number of time periods,

Xj = {Xjl}' t = 1, ... ,T is a variable that enters a criterion of a type SUP, INF, FOL or DER.

Transformation will be discussed for each type of criteria:

Type: MIN

additional row (with name which is concatenation of following three strings:

, <

',eni, , ... ' ) is generated in form:

-v

+

wiqi

S

wiiii

Type: MAX

is transformed in the way similar to type MIN, with additional (internal, for computations only) change of the signs ofWi to negative.

Type: SUP

additionalT rows (with names which are concatenations of strings'

<

',eni, '.'eht, where t = 1, ... ,T) are generated in forms:

-v

+

WiXjl ~WiXjl

+

Wiiji Type: INF

is transformed in the way similar to type SUP, with additional (internal, for com- putations only) change of the signs ofWi to negative.

Type: FOL

• additional T columns (with names which are concatenations of strings'

+ "

eni, '.', chI,wheret = 1, ...,T) are generated; in the following formulae this name is replaced by e~

• additional T columns (with names which are concatenations of strings' - "

eni, '.', chI,where t

=

1, ...,T) are generated; in the following formulae this name is replaced by ei;

• additionalT rows (with names which are concatenation of strings'

= "

eni, '.', chI, wheret

=

1, ... ,T) are generated in form:

e~ - ei; - Xjl

=

^-Xjl

• additionalT rows (with names which are concatenations of strings'

<

',eni, ,.', chI, wheret = 1, ...,T) are generated in the form:

-v

+

wi(e~

+

^ei;) ~ ^Wiiji Type: DER

• additional 2X T columns are generated in the same way as described for a criterion of the type FOLj

• additionalT rows (with names with are concatenations of strings'= " eni, '.', chI, where t = 1, ...,T) are generated in form:

e~

- e;; -

^Xj,l

+

Xj,l-l = O.

• additional T rows (with names which are concatenations of strings'

<

',eni, '.', chI) are generated in form:

-v

+

wi(e~

+

ei;) ~ ^Wiijj

Auxiliary goal function, which is to be minimized, is generated in the following form:

v

+ _~m (2.:

^Wiqi

⁺ 2.:(2.:

^WjXjl

⁺ 2.:

^WI;(etl

⁺

ekl )))

i t ^j I;

where summation is done over corresponding sets of respective criteria, Le. indices i,j,k correspond to criteria of type: MIN or MAX, SUP or INF and FOL or DER, respectively; ~mis given parameter.

The name of auxiliary variable v is '.. dummy.', whereas the name of auxiliary goal function is '.dummy..'.

Value of~m may be changed by the command MEPS in a routine for modification of multicriteria parameters.

2.5 Formulation of single criteria linear-quadratic problems

HYBRID 3.1 allows for solution of a single-criterion linear-quadratic problem with a simple quadratic term. For a problem which does not have recognized structure (as discussed in sect. 2.3) the formulation takes the following form:

minex

+ b/2)llx - xl1

²

subject (2) and (3), where

x

is a given point in the solution space and 1

>

0 is a given parameter.

Similarly, for a dynamic problem one may formulate the problem in the following way:

T

min Latxt

+

bu

+ b/2)llu - ul1

²

t=]

subject (9) and (11), where

u

is a given point in the space of independent variables.

3 Theoretical foundations and problems

3.1 General remarks

The most popular methods for solving linear programming problems are based on the simplex algorithm. However, a number of other iterative non-simplex approaches have recently been developed (Mangasarian, 1981, Polyak and Tretiyakov, 1972, Sosnowski, 1981). HYBRID belongs to this group of non-simplex methods. The solution technique is based on the minimization of an augmented Lagrangian penalty function using a modification of the conjugate gradient method. The Lagrange multipliers are updated using a modified version of the multiplier method (Bertsekas, 1976) (see Sections 4.2 and 4.4).

This method is useful not only for linear programming problems but also for other purposes, as described in Section 1.2. In addition, the method may be used to solve problems with non-unique solutions (as a result of regularization-see Section 4.5).

The following notation will be used:

ai - denotes the i-th row of matrix A

Xi - denotes the j-th component of vector x

II xii -

denotes the Euclidian norm of vector x

(u)+ - denotes the vector composed of the non-negative elements of vector u (where negative elements are replaced by zeros)

AT - denotes the transposition of matrix A.

3.2 The multiplier method

We shall first explain how the multiplier method may be applied directly to LP problems.

Consider the problem (PO), which is equivalent to the problem (P) defined in Sec- tion 2.1:

min ex

Bx~ d (PO)

where d E RP, B is a p X n matrix, and m ~p _~ 2(m

+

n). To apply the multiplier method to this problem we proceed as follows:

Select initial multipliers yO (e.g., yO = 0) and pER, p

>

O. Then for k = 0,1, ...

determine successive values ofXk+l, yk+l where

and where

L(x,ylc)

=

ex

+ (II

(ylc

+

p(Bx - d))+112

-

Ily

^1c112) / (2p) until a stopping criterion is satisfied.

The method has the following basic properties:

(13) (14)

1. Apiecewise quadratic differentiable convex function is minimized at each iteration.

2. The algorithm terminates in a finite number of iterations for any positivep.

3. There exists a constant15such that for any p ~15the algorithm terminates in the second iteration.

Note that it is assumed above that the function L(·,ylc) is minimized exactly and that the value of the penalty parameterpis fixed. Less accurate minimization may be performed provided that certain conditions are fulfilled (see, e.g., Sosnowski, 1981, Bertsekas, 1976). For numerical reasons, a non-decreasing sequence of penalty parameters {pic}

is generally used instead of a fixed p.

3.3 The conjugate gradient method for the minimization of an augmented Lagrangian penalty function

The augmented Lagrangian function for a given vector of multipliers y will be called the augmented Lagrangian penalty function (Flecher, 1981). For minimization of that function the conjugate gradient method has been modified to take advantage of the formulation of the problem. The method may be understood as an modification of the techniques developed by Polyak (1969), O'Leary (1980) and Hestenes (1980) for minimization of a quadratic function on an interval using the conjugate gradient method.

The problem (P) may be reformulated as follows:

mmex Ax

+

z = b

(PS)

where z E Rm are slack variables.

Formulation (PS) has a number of advantages over the initial formulation (PO):

1. The dimension of matrixA in (PS) is usually much smaller than that of matrixB in (PO).

2. The augmented Lagrangian problem is one of minimization of a quadratic function in (PS), and of minimization of a piecewise quadratic in (PO).

3. Some computations only have to be performed for subsets of variables. Note that slack variables are introduced only for ease of interpretation and do not have to be computed.

In (PS) the augmented Lagrangian is defined by

L(x,z,U) = ex

+ (Ilu +

p(Ax

+

z -

b)11 2-lluI1 2) I

(2p). (15) We shall first discuss the problem of minimizingL(x, z, U) for givenU,p

>

0,subject to lower and upper bounds for x and z. Let us consider the following augmented Lagrangian penalty function

F(x, z) = (el p)x

+ (Ilul

p

+

Ax - b

+ zl12 - Ilulp112)

12.

The gradient ofF is defined by

aFax =elp+AT(z-g)

- = z - gaF az

where

(16)

g

= - u I

p - Ax

+

b.

From the Kuhn-Tucker optimality condition, the following relations hold for the minimum point (x·, z·):

ar if xi = I;, ar

if xi = u;,

- > 0 - < 0

ax; - ax; -

ar if aF·

if

- > 0 z;" = 0, - < 0

z; ==

^Ti,

aZ_i - aZ_i -

and

aF·

- - =

0 if I;

<

xi

<

u;

ax;

aF·- - = 0 if 0

<

z;"

<

ri.

aZi

For any given point such that I ~ x ~ ^U it is possible to determine slack variables

o

~ Z ~ r in such a way that the optimality conditions with respect to z are obeyed.

Variables z are defined by

Zi =

1 ~i

gi

if gi _~ 0 (aFjaz;

>

0) if gi;::: ri (aFjaz;

<

0) if ri

>

^gi

>

0 (aFjaz; = 0).

(17)

AI AI

J J

Al Al

J J

We shall use the following notation and definitions. The vector of variables x with indices that belong to a set J will be denoted by x^J, and analogous notation will be used for variablesg. Letqdenote minus the gradient of the Lagrangian penalty function reduced to x-space(q = - (aF j ax)). The following sets of indices are defined for a given point x:

The set of indicesI of violated constraints, i.e.,

1= { i :gi ;::: ri}U { i :gi _~ O}.

1

is the complement ofI, i.e.,

1=

{1,2, ... ,rn}\I.

The set of indices I can be also interpreted as a set of active simple constraints for z.

The set of indices J of variables that should be equal to either the upper or the lower bound, Le.,

J

=

{j :^Xj

=

Ij and qj ~ O} U{j :^Xj = Uj and qj;::: O}.

J

is the complement ofJ, i.e.,

J =

{l,2, ...,n}\ J.

For the sake of illustration the matrixA may be schematically split up in the follow- ing three ways (see the Figure below): first according to active rows, second according to basic columns and third with illustrate the part of the matrixA for which augmented Lagrangian penalty function is computed. The contents of the matrixA~(for which the augmented Lagrangian penalty function is computed) changes along with computations.

J

J

EE

In essence, the augmented Lagrangian penalty function is minimized using the con- jugate gradient method with the following modifications:

1. During the minimization process x and z satisfy simple constraints and z enters the augmented Lagrangian in the form defined by (17).

2. The conjugate gradient routine is run until no new constraint becomes active, Le., neither set I nor set J increases in size. Ifthis occurs, the computed step length is shortened to reach the next constraint, the corresponding set (lor J) is enlarged and the conjugate gradient routine is re-entered with the direction set equal to minus the gradient.

3. Sets J and I are defined before entering the procedure discussed in point 2 and may be only enlarged before the minimum is found. When the minimum with respect to the variables with indices in sets

J

and I has been found, sets J and I are redefined.

4. Minimization is performed subject only to those components of variables x whose indices belong to set

J,

i.e., variables that are not currently equal to a bound value.

5. Minimization is performed subject only to those components of variables z whose indices do not belong to set I, i.e., slack variables that correspond to non-active simple constraints for z. Note that, formally, this requires only the use of different formulae for z. In actual fact it is sufficient to know only the set I, which defines the minimized quadratic function.

4 Solution technique

4.1 Algorithm for minimization of augmented Lagrangian

We may now present the algorithm for minimization of the augmented Lagrangian penalty function in a more formal way. The algorithm consists of the following steps:

1. For given y and p

> a

choose a point x such that l ~ x ~ u 2. Computeg

=

-.Y/p - Ax

+

^b

3. Determine sets I and

I

1= {i :gi

>

Ti } U {i :gi

<

a}, 1={l, ... ,m}\I

4. DefineYas follows:

gi - Ti

Yi

^{= {}

gi

5. Compute the minus gradient:

if gi - Ti

> a

otherwise

6. Determine sets J and

J

J

=

{j :^Xj

=

Ij and qj :S o}U{j :^Xj

=

^Uj and qj ~ 0 }

J={I, ... ,n}\J

7. Ifqj

=

0 for all j E

J

then x is a minimum point of the augmented Lagrangian penalty function

8. Set pJ = qJ 9. Compute

S

=

^AJpJ

h=

IIll1

²

d=

IIs

^I

1I

² a(l)

=

hid

Note that a(l) is the conjugate gradient step length in direction pJ

10. Find the step length that would violate the nearest non-active constraint, Le., for i E1,

a(2)

=

~W{gds;}, K

=

{i: i E1,si

>

o}

a(3) = ~W{(gi - Ti)ls;}, K = {i: iE 1,si

<

o}

11. Find the step length that would enable a variable to reach a bound, i.e., a(4) = min_jEK (I· - x·)lp·

J J J'

a(5)

=

min(u· -_jEK x·)lp·

J J J'

K

=

{j :j E J,Pj

<

o}

K = {j :j E

J,

Pj

>

0 }

12. Determine step length

a =

min;=l, ... ,s(a(i)). If

a =

min(a(2),a(3)) add the row index for which this condition holds to set I and remove that index from set

1.

If

a =

min(a(4),a(5)) add the column index for which this condition holds to set J and remove that index from set

J.

13. Compute the new point xJ :=x^J

+

apJ and the minus gradient at that point:

gi :=gi -

aS

i

14. IfqJ = O. go to step 2

15. If

a =

.a(l) continue with the conjugate gradient step, i.e.

{J =

IIq

^J

I1

²

I

h pJ :=

l +

{JpJ and go to step 9

16. Go to step 8

Note that the condition qJ gradient tolerance.

o

is in practice replaced by

IlqJ11

~ ^E:, where E: is a

4.2 Steps of the multiplier method

Let the violation of i-th constraint in a point x^k be defined in the following way:

vt

=

max{aixk - bi ,-a;xk

+

bi - Ti'O}

and

Ilvkll

oo denotes the100 norm of violated constraints. The multiplier method will be presented in algorithmic form.

1. Compute an initial vector of multipliers on the basis of the particular option chosen (Le., either yO

=

0 or yO corresponding to the constraints violated at starting point x)

2. Find x

Hl

which minimizes the augmented Lagrangian penalty function (see Sec- tion 3.3) with accuracy E:^k• It is assumed that

where the sequence E:k

---+

o.

In addition, E:mi :::: E:k

:::: E: m:>:, where E:mi. E: m:>: is the assumed minimum and maximum accuracy, respectively.

3. Compute new multipliers

if yk

+

pk(ai xHl - b;) :::: 0 if yk

+

l(a;x^{Hl -} b;

+

^T;) ~0 otherwise

4. If

IlyHl-yk II >

^E:dthen setpk+1

=

min(pk P"Pm:>:) , P.

>

1., Pm:>:is a given maximal value of the penalty parameter.

Set E:H 1= E:kE: •• where E:.

<

1. is an assumed parameter.

Set k ^:= k

+

1 and go to step 2

5. Set k := k

+

1 and find ^XH1 which minimizes the augmented Lagrangian. IfXH1

is feasible (

Ilvkll

~ F EAS ) then assume it as a solution and stop.

Otherwise setpHI

=

min(pk p" Pm:>:), and E:k+1

=

^E:kE:. and go to step 2.

4.3 Solution technique for DLP

We will not repeat reasoning given in the first part of sect. 2.3. Instead, let us point out basic differences between the algorithms for static LP and DLP:

1. Minimization is reduced to a subspace of decision variables. Gradient of La- grangian penalty function is computed for variables that belong to a subspace of decision variables. This (together with arguments already presented in sect.

3.3) shows advantages due to the use of dynamic structure of DLP problem in comparison with presentation of such a problem as a large LP.

2. The structure of matrices BI , ••• ,BT and Fa, ... , FT has no impact for the algo- rithm nor affects the technique of storage of data, because super-sparse technique is applied (cfsect. 1.4). It should be also pointed out that the method of trans- forming a multicriteria problem to a parametric LP one introduces constraints (cf sect. 2.4.3) that-for the proposed (cfsect. 2.4.2) types of criteria-do not fit to the staircase structure of CDLP (cfPropoi, 1976). Therefore, any technique that would exploit the staircase structure of DLP would also imply a reduction of a number of criteria types. The alternative is then to treat a problem as a large LP static one or to apply a technique that does not exploit the classical DLP structure.

3. State equations are solved (for given decision variablesu) by forward substitution.

Therefore any single constraints for state variables have to be treated as general constraints and included into the matrix. Gradient need not to be computed for those variables, but state equation is solved twice (for state variables and variations) .

4. A conjugate trajectory 'l1 is computed from conjugate equation by backward sub- stitution and has an interpretation of dual variables for state equations. No other variables associated with those rows (defined in sect. 3.3, i.e. Lagrange multipliers, shifted constraintsg) are computed for state equations rows.

5. The general structure of the algorithm for DLP is similar to that presented in sect. 3.4. To sum up basic differences one may observe that:

we consider a problem that is equivalent to a static LP but reduced to the subspace of decision variables and is solved in the way similar to that de- scribed in sect. 3.3 and 3.4,

state equations are solved for control variables and for variations, a conjugate trajectory 'l1 is computed.

4.4 Algorithm for minimization of augmented Lagrangian for DLP

Now we may present the algorithm for minimization of the augmented Lagrangian function for DLP in a more formal way. In each iteration of multiplier method, the

following optimization problem is solved: minimize the augmented Lagrangian penalty function

T

F(x, u, z)

= 2:(al/

^P)XI

+

(b/p)u+

1=1

subject to

I-I

-Hlxl

+ 2:

^AI-I,iXi

+

Blu= CI

;=0

with a given initial condition Xo and

t = 1, ...,T

O~Z~T

where z is a vector of slack variables, which-as discussed in sect. 3.3-are not used in the algorithm. The algorithm consists of the following steps:

1. For given y and p choose a point u such that e~ u ~

f

2. Solve the state equation

I-I

Hlxl

= 2:

^AI-I,iXi

+

Blu - CI i=O

with given initial condition Xo

3. Compute shifted constraints for constraints (10)

t

=

1, ... ,T

T

g

= -y/

^{p -}

2:

^{Flxl - Du}

+

^d

1=0

and determine sets I, j

I

=

{i :gj

>

^Ti}U { i :^gi

<

0 } while j is the complement ofI.

gi - Tj

4. Definegas follows :

gi

= {

gi

if gi

>

^Ti

otherwise

5. Find the conjugate trajectory by solving backwards the conjugate equations

HTw

I = T-I

2:

ATtWt+!

+

(F/)TgI - at/p,

i=1

with boundary condition

t = T - l , ... ,1

6. Compute the minus gradient reduced to subspace of decision variables

T q= -b/p

+

(D¹

f

^g1

+ L

^B'{Wt

t=l

7. Determine sets J and

J

J

=

{j :u;

=

f; and q; ~ 0 } U{j :u;

= /;

and q; ~ 0 } while

J

is the complement ofJ

8. Ifq; = 0 for all j E

J

then u is a minimum point of the augmented Lagrangian penalty function

10. Solve state equation in variations

t = 1, ... ,T

with boundary conditionao

=

0 11. Compute

T S = D^JpJ

+ L

^Fiat

t=o

h=

IIll1

²

V

= IIs

¹

11

²

a(l)

=

h/v

Note that a(l) is the conjugate gradient step length in directionpJ

12. Find the step length that would violate the nearest non-violated constraint, i.e., a(2) = ~W{gds;}, K = {i : i ElandSi

>

o}

a(3)

=

rp.iKn{(gi - Ti)/S;}, K

=

{i: i Eland Si

<

o}

.E

13. Find the step length that would enable a variable to reach a bound, i.e.,

a(5)

=

~W{(/;

-

u;)/p;},

K = {j : j E

J

andp;

<

0 } K = {j : j E

J

and p;

>

o}

14. Determine step length

a

^{= .min (0: (i) )}

1=1 •...,5

If

a ⁼

min( 0:(2),0:(3)) add the row index for which this condition holds to set I and remove that index from set

1.

If

a

⁼ min(0:(4), 0:(5)) add the column index for which this condition holds to set J and remove that index from set

J.

15. Compute:

u/:= u^J

+

apJ

Xl :=Xl

+

aO"I

gi := g; - as;

16. For the new gl solve the conjugate equation (as in step 5) 17. Compute the minus gradient:

l =

_bJjp+ (D~fgl

+ L(BdjWl

T

1=1

18. If qJ

=

0, then go to 2

19. If

a =

0:(1) continue with the conjugate gradient step, i.e.

{3

⁼

IIll1

²

j

h

and go to step 10 20. Go to step 9

Note that the condition qJ = 0 is in practice replaced by

IlqJ11

~ c. The value of c may be quite large in the first few iterations; it then decreases as the number of iterations increases.

4.5 Regularization

Itis possible that a linear programming problem may have nonunique optimal solutions.

Although this is theoretically rare, in practice many problems actually have a large set of widely varying basic solutions for which the objective values differ very little (Sosnowski, 1981). In some cases, the simplex algorithm will stop when a basic solution is recognized as optimal for a given set of tolerances. For problems with a nonunique optimum, the first optimal solution found is accepted, so that one may not even be aware of the non-uniqueness of the solution reported as optimal.

Thus we are faced with the problem of choosing an optimal (or, in most cases, to be more accurate, a suboptimal) solution that possesses certain additional properties required by the user. This problem may be overcome by applying an approach called