Graphical Representation of
Constrained Optimization Problems

The author describes a function for automatically generating plots to demonstrate optimization solutions.

by M. Asghar Bhatti

Constrained optimization problems arise naturally in many different disciplines. A structural engineer designing a multistory building must choose materials and proportions for different structural components to create a safe structure that is as economical as possible. A portfolio manager for a large mutual fund company must choose investments that generate the largest possible rate of return for its investors while keeping the risk of major losses to acceptably low levels. A plant manager in a manufacturing facility must schedule the plant operations such that the plant produces products that maximize company's revenues while staying within the available resource limitations. All these situations can be formulated in terms of a constrained optimization problem that is expressed in the following form.

Find a vector of optimization variables, x=((((Subscript[x, 1]), (Subscript[x, 2]),  ... , (Subscript[x, n])))^T)in order to

Minimize an objective function, f(x)

Subject to

({{(Subscript[g, i])(x)≤0, \(\(i\  = \ 
                1\)\(,\)\(\ \)\(2\)\(,\)\(...\)\(,\)\(\ \)\(m\)\(	\)\), Inequality constraints}, {(Subscript[h, i])(x)=0, i  =  
                1,  2,   ... ,  p, Equality constraints}})

Many sophisticated methods exist for solution of these problems and are described in a variety of books on the subjec; for example see [Bhatti 1999]. For problems involving two variables, a graphical representation is frequently used in explaining these solution methods. A typical graph is shown in Figure . The dark circle represents the equality constraint (h((Subscript[x, 1]), (Subscript[x, 2])) = 0). The straight line represents the boundary of the linear inequality constraint (g((Subscript[x, 1]), (Subscript[x, 2])) = 0). On one side of this line g((Subscript[x, 1]), (Subscript[x, 2])) < 0 and thus any point on that side satisfies the inequality constraint. On the other side g((Subscript[x, 1]), (Subscript[x, 2])) > 0 and therefore represents the infeasible side of the constraint. On the graph the infeasible side of an inequality constraint is shown by lightly shading that side.

A feasible domain is the region over which all constraints are satisfied. Considering both constraints together, we see that the feasible domain for the example graph is the part of the circle h that is below the line g.  Our next task is to identify one or more feasible points where the objective function is minimized. Contours of the objective function help in making this determination. The thin gray lines, labelled -2, -0.66, and 2, represent objective function contour lines. The curve with f = -2 lies entirely outside the feasible domain. The curve with f = 2 intersects the feasible domain in three points. It is clear that a point that is feasible and also has the lowest value of the objective function is the one where the objective contour just touches the feasible domain. In the figure this occurs for f = -0.66, and this the optimum solution.

The next section describes the basic steps necessary for generating the figure shown. The steps are implemented in a function called GraphicalSolution. This function and its various options are shown in the following section. The use of some of these options is illustrated through several examples. The main ideas used in the implementation are simpl,e but the actual Mathematica code is fairly lengthy, primarily because of all different options that have been provided to give a satisfactory graph for a wide variety of optimization problems. Therefore the actual code is not presented in the paper. The complete code is available in the electronic supplement.

Basic Steps

The graph shown in Figure  represents the solution of the following optimization problem.

Minimize	f  =   (Subscript[(e^x), 1]) -  (Subscript[x, 1]) (Subscript[x, 2])  +  (Subscript[x, 2]^2)

Subject to	(({{g  =  2  (Subscript[x, 1])  +  (Subscript[x, 2])  -  2 ≤ 0}, {h  =  (Subscript[x, 1]^2)  +  (Subscript[x, 2]^2)  -  4 = 0}}))

As most experienced Mathematica users know, all graphics functions do not gracefully accept subscripted variables. To overcome this difficulty, we define the problem in terms of non-subscripted variables as follows.

f = (E^x1) - x1 x2 + (x2^2);
g = 2 x1 + x2 - 2;
h = (x1^2) + (x2^2) - 4;

The ContourPlot function is the main tool that is used to create the graphical representation. The constraint boundaries are shown by drawing contours with function values equal to zero. For the objective function we simply need to specify the number of contours: The ContourPlot function automatically determines the appropriate contour values. Thus the graph shown in Figure  can easily be obtained as follows.

Block[{$DisplayFunction = Identity}, 
        hp = ContourPlot[h, {x1, -4, 4}, {x2, -4, 4}, 
            Contours -> {0}, ContourShading -> False, 
            ContourStyle -> {{Thickness[0.01]}}, 
            PlotPoints -> 50]; 
gp = 
          ContourPlot[g, {x1, -4, 4}, {x2, -4, 4}, 
            Contours -> {0}, 
            ContourShading -> False]; 
fp = 
          ContourPlot[f, {x1, -4, 4}, {x2, -4, 4}, 
            Contours -> 10, ContourShading -> False, 
            PlotPoints -> 50, 
            ContourStyle -> {{GrayLevel[
                    0.7]}}];
];
\(\(Show[{fp, gp, hp}, AspectRatio -> Automatic, Axes -> True, 
        AxesLabel -> {\*""(Subscript[x, 1])"", \*""(Subscript[x, 2])""}];\)\)

This graph has three main shortcomings.

1. The constraints are not labelled. For this simple example, it is not a major problem because there are only two constraints. The equality constraint is shown with a thicker line and thus is easy to identify. However if there were several constraints, it would be necessary to identify constraints by placing appropriate labels next to their contours.

2. The objective function contours are not labelled and therefore there is no way to identify the value represented by a particular contour. Because of this limitation it is impossible to determine the optimum point from the graph without performing additional computations.

3. There is no indication of the feasible side of the inequality constraint.

In order to place labels for constraint functions we can use the Text function and place suitable text at the appropriate locations. Manually the text locations can be determined easily from the graph. However automatically determining label locations requires knowledge of lines used to draw the contours. One way to accomplish this is to first convert the ContourGraphics object returned by ContourPlot to a graphics object by wrapping it with Graphics. Now the data related to Lines can be selected using the Cases function. From the arguments of Line function one can then select possible point for text placement. As an example, for the inequality constraint, the first point of contour line is selected as follows.

gpLoc = Flatten[
          Cases[Flatten[Graphics[gp][[1]]], 
            Line[pts_] -> pts]][[{1, 2}]]

The same procedure can be used for all other constraints in a problem. A similar procedure can also be used to determine objective function contour values and then these values can be shown at appropriate locations. However it is not necessary because an efficient implementation for labelling contours is widely available in the ExtendGraphics package by [Tom Wickham-Jones, 1994]. The package is called LabelContour and the function is called LabelContourLines.

(*Glenn. The name Ruskeepaa below is supposed to end with unlaut a [twice]. I can't get that character to stick in this notebook. Something is weird. David*)

Automatic shading of the infeasible side is the last hurdle left. For linear constraints it is not very difficult to determine the feasible side, and perhaps can be shown by polygon shading; for example, see [Heikki Ruskeepää 1999]. However, for nonlinear constraints the implementation presented in this paper is the only one that provides an automatic solution. The main idea is to draw a second contour representing g((Subscript[x, 1]), (Subscript[x, 2])) = ϵ where ϵ is a small positive number. This contour then will be on the infeasible side of the constraint boundary. The infeasible contour is drawn with a different thickness and shading. By adjusting the ϵ and thickness one can obtain graphs that automatically represent infeasible side in a satisfactory manner.

Combining all these steps, as follows, we obtain the graph shown in Figure .

Needs["ExtendGraphics`LabelContour`"];

\(\(Block[{$DisplayFunction = Identity}, 
          hp = ContourPlot[h, {x1, \(-4\), 4}, {x2, \(-4\), 4}, 
              Contours -> {0}, ContourShading -> False, 
              ContourStyle -> {{Thickness[ .01]}}, 
              PlotPoints -> 50]; 
hpLoc = \(Flatten[
                Cases[Flatten[\(Graphics[hp]\)[\([1]\)]], 
                  Line[pts_] -> pts]]\)[\([{1, 
                  2}]\)]; 
gp = 
            ContourPlot[g, {x1, \(-4\), 4}, {x2, \(-4\), 4}, 
              Contours -> {0}, 
              ContourShading -> False]; 
gpe = 
            ContourPlot[g, {x1, \(-4\), 4}, {x2, \(-4\), 4}, 
              Contours -> {0.38}, ContourShading -> False, 
              ContourStyle -> {{GrayLevel[ .8], 
                    Thickness[ .03]}}]; 
gpLoc = \(Flatten[
                Cases[Flatten[\(Graphics[gp]\)[\([1]\)]], 
                  Line[pts_] -> pts]]\)[\([{1, 
                  2}]\)]; 
fp = 
            LabelContourLines[
              ContourPlot[f, {x1, \(-4\), 4}, {x2, \(-4\), 4}, 
                Contours -> 7, ContourShading -> False, 
                PlotPoints -> 50, 
                ContourStyle -> {{GrayLevel[ .7]}}], 
              LabelPlacement -> Automatic];
];\)

      \(Show[{fp, gp, gpe, hp}, AspectRatio -> Automatic, 
          Axes -> True, 
          AxesLabel -> {\*""\!\(x\_1\)"", \*""\!\(x\_2\)""}, 
          TextStyle -> {FontFamily -> "Times", 
              FontSize -> 9}, 
          Epilog -> {Text["h", hpLoc + { .1,  .15}], 
              Text["g", gpLoc + { .15, 0}]}];\)\)

GraphicalSolution function

The ideas presented in the last section form the basis for a function called GraphicalSolution that is designed to automatically generate graphical representation of optimization problems. The function is implemented to handle both constrained and unconstrained problems. Several options are implemented to effectively deal with a variety of problems. The process of placement of labels for constraints and identification of infeasible side of inequality constraints is automated. The only manual activity is the identification of the optimum point, which obviously is impossible to automate without using a sophisticated constrained optimization algorithm. The actual placement of the dark circle representing the optimum point and the word "Optimum" is accomplished by using the standard Epilog option available in the Show function.

Needs["OptimizationToolbox`GraphicalSolution`"];

?GraphicalSolution

The function allows subscripts on the variable names. Since all built-in graphics functions do not accept subscripted variables, internally the optimization variables are mapped to local variables without subscripts. The original variables are used in the axes labels, however.

The GraphicalSolution function accepts the following options.

OptionsUsage[GraphicalSolution]

The last three options provide further useful educational tools. The gradient vectors at selected points can be drawn to show physical interpretation of optimality conditions for constrained problems. The PlotHistory  option superimposes a line joining the given points on the graph. This option is useful to show the progress of numerical algorithms towards optimum.

In addition to the above options, most relevant options for ContourGraphics and Graphics functions are passed on to these functions. Epilog option in particular is very useful and is used in the following examples.

Example 1

This example shows a problem involving an equality constraint. Only the points lying on the dark line representing constraint h are feasible.

Minimize	f  =   (x^2)  +  (y^2)
Subject to	h  =  
        3 (x^2)  +  4 x y  +  6 (y^2)  -  140

f = (x^2) + (y^2);
h = 3 (x^2) + 4 x y + 6 (y^2) == 140;

grph = 
        GraphicalSolution[f, {x, -20, 20}, {y, -20, 20}, 
          Constraints -> {h}, 
          ObjectiveContours -> Table[i, {i, 20, 80, 20}], 
          PlotPoints -> 50];

As seen from this graph we have following two global minimum points.

xopt1 = {2, 4}; xopt2 = {-2, -4};
f /. {ToRules[{x, y} == xopt1], 
        ToRules[{x, y} == xopt2]}

Show[grph, 
        Epilog -> {RGBColor[1, 0, 0], Disk[xopt1,  .2], 
            Text["  Optimum", xopt1, {-1, 0}], Disk[xopt2,  .2], 
            Text["  Optimum", xopt2, {-1, 0}]}];

Example 2

This example demonstrates the options useful for geometric interpretation of Karus-Kuhn-Tucker (KT) optimality conditions. These conditions essentially state that, at a minimum point, the gradient of the objective function is numerically equal to a linear combination of constraint gradients but points in the opposite direction, see [Bhatti, 1999] for details.

f = (x^2) - 4 y + (y^2);
g = {-(x^2) + (y^2) ≤ 0, 
          1/4 - Log[x + 2  y] ≤ 0};
vars = {x, y};

Using the GradientVectors option, the gradients of f and (Subscript[g, 1]) are shown at point (1, 1). The graph clearly demonstrates that this point satisfies KT conditions.

\(\(GraphicalSolution[f, {x,  .1, 2}, {y,  .1, 2}, Constraints -> g, 
        ObjectiveContours -> 5, ShadingOffset -> 0.05, 
        ShadingThickness ->  .015, 
        GradientVectors -> {{f, 
              First[g[\([1]\)]]}, {"∇f", \*""∇(Subscript[g, 1])""}, {{1, 1}}}, GradientVectorScale -> 0.1, 
        Epilog -> {RGBColor[1, 0, 0], Disk[{1, 1},  .03], 
            Text["Minimum", {1.1, 1}, {\(-1\), 0}]}];\)\)

Using the GradientVectors option once again, the following graph demonstrates that the intersection of constraints (Subscript[g, 1]) and (Subscript[g, 2]) does not satisfy the KT conditions.

sol = NSolve[Thread[Map[First, g] == 0], vars]

\(\(GraphicalSolution[f, {x,  .1, 1}, {y,  .1, 1}, Constraints -> g, 
        ObjectiveContours -> 5, ShadingOffset -> 0.018, 
        GradientVectors -> {{f, First[g[\([1]\)]], 
              First[g[\([2]\)]]}, {"∇f", \*""∇(Subscript[g, 1])"", \*""∇(Subscript[g, 2])""}, {{x, y} /. 
                sol[\([2]\)]}}, GradientVectorScale -> 0.1];\)\)

Example 3

This example features the PlotHistory option that is used to superimpose a solution path taken by a numerical optimization algorithm on the graph. The problem is defined as follows.

f = Exp[-(Subscript[x, 1]) - (Subscript[x, 2])] + (Subscript[x, 1]) (Subscript[x, 2]) + (Subscript[x, 2]^2);
g = {1/4 Exp[-(Subscript[x, 1])] + 1/4 (Subscript[x, 2]^2) - 1 ≤ 0, 
          (Subscript[x, 1])/5 + (Subscript[x, 2])/2 - 1 ≤ 0};

The problem has been solved using a constrained optimization algorithm that generated the following history of points.

xhist = {{1, 1}, {0.7838, 0.2828}, {0.7988, 
            0.03191}, {1.203, -0.395}, {2.044, -0.3617}, {2.591, -1.496}, {3.506, -1.128}, {4.117, -1.707}, {5.915, -2.029}, {7.966, -2.000}, {9.970, -2.000}, {10.00, -2.000}};

\(\(GraphicalSolution[f, 
{x\_1, \ \(-3\), 12}, \ {x\_2, \ \(-3\), 5}, 

        ObjectiveContours -> 5, Constraints -> g, \ ShadingOffset -> 0.045, 
        PlotPoints -> 35, 
        Epilog -> {RGBColor[1, 0, 0], Disk[Last[xhist], 0.1], 
            Text["Minimum", Last[xhist] + {0, \(- .2\)}], 
            Text["Start", First[xhist] + {0,  .2}], 
            Line[xhist]}];\)\)

Concluding  remarks

A function called GraphicalSolution is presented for automatically generating plots useful for demonstrating optimization concepts and solution algorithms. The basic idea of the function is explained. Several examples are presented to demonstrate many options implemented for the function.

References

ABOUT THE AUTHOR

Asghar Bhatti received his B.E. degree in Civil Engineering from the University of Karachi, Karachi, Pakistan, in 1973. After working in Pakistan for a year in a large consulting engineering firm, he attended the University of California, Berkeley, where he received his M.S. degree in 1975 and Ph.D. in 1979, both in Structural Engineering and Structural Mechanics. He has published several papers and reports in this area, and is the author of  Practical Optimization Methods with Mathematica Applications.  New York: TELOS, Springer-Verlag, 1999. His teaching interests include finite element techniques in engineering, advanced structural analysis, analysis and design of structures subjected to dynamic loads, optimal design, application of computers in structural engineering and computer graphics.

Asghar Bhatti
Department of Civil and Environmental Engineering
University of Iowa
Iowa City, IA 52242
mabhatti@uiowa.edu

ELECTRONIC SUBSCRIPTIONS

Included in the distribution for each electronic subscription is the file bhatti99.nb, containing Mathematica code for the material described in this article.