# The MATLAB Notebook v1.5.2

• Doc File169.00KByte

• ### Solving a system of equations in matlab

• Doc File 169.00KByte

﻿

Section 7: More on the Logarithmic Barrier Method in Linear Programming.

In the previous section, we introduced the logarithmic barrier method as an example of an interior method in linear programming. As it has been presented so far, this method depends on the ability of MATLAB to find the critical points of the modified objective function by symbolically solving equations. With n variables and k constraints, this amounts to solving a system of n polynomial equations of degree k, which rapidly exhausts the capacity of MATLAB.

In the present section we present an alternative procedure, which is essentially a higher dimensional form of Newton's method for finding the roots of an equation. We begin by reviewing the problem of finding roots of equations involving one variable. To solve the equation ax + b=0, one simply sets [pic]. Obvious as this is, it is the basis of Newton's method. For any differentiable function f, we have, to first order in x-x0, f(x)=f(x0)+ (x-x0)f ' (x0), which gives for f(x)=0, the approximate solution [pic]. The closer x0 is to an actual root of f(x)=0, the better the approximation is. Newton's method is simply to iterate this idea, with x1 taking over the role of x0.[pic]Let us see how this works to find a root of the equation ex - 3 -x =0. We will start with x0 = 1, and see what happens.

syms x

f=exp(x)-3-x

f =

exp(x)-3-x

x0=1

x0 =

1

This initializes the process. Now each iteration of the next input block will bring us closer to an actual solution. The second line measures how close we are, so we can decide when we are close enough.

x1=x0-subs(f,x,x0)/subs(diff(f),x,x0)

subs(f,x,x1)

x0=x1

This tends to work well as long as we start reasonable near an actual solution and are not unlucky enough to encounter a point at which f '=0, which would cause the first line to fail, or f' very small, which would tend to cause serious overshoot. It is for this reason that we began with x0=1 instead of x0=0. Of course we can find points at which f '(0) = 0 in much the same manner, applying Newton's method to the function f '.

When we come to functions of several variables, we do not expect isolated solutions to equations like f(x,y,z)=0. In order to have generically isolated solutions, we need the same number of equations as variables. An important special case of this, and the one we are concerned with in connection with the logarithmic barrier method, is finding critical points for a function of several variables. Since there are as many partial derivatives as there are variables, critical points are generically isolated.

To find a critical point of a function of several variables by the analogue of Newton's method, we need to start at a point P0, and look at f(P) near P0 to second order in P-P0. Recall from Section 4 of these notes (with a slight change of notation) that, to second order in P-P0, we have

[pic] Here, we recall that Hf denotes the Hessian matrix of f, and [pic] denotes the gradient of f. We are writing the coordinates of our points horizontally. Since the components of the Hessian matrix are the partial derivatives of the components of the gradient, we also have, to first order in P-P0, [pic]. If Hf=H were constant, the solution of [pic] would be [pic], as may be verified by substitution. This is the basis for our adaptation of Newton's method. In the general case P1 will not be the desired solution but we hope it will be an improvement on P0, and we can iterate as before.

Unfortunately, if we try to apply this idea directly in the logarithmic barrier method, we will probably not only overshoot the desired solution, but jump out of the feasible region altogether. For this reason, we modify the method by using [pic] to determine the direction in which we should move, but directly specifying the distance we will move, by adding to P0 a controlled multiple of a unit vector in the desired direction. At every stage, we check to make sure we are still in the feasible region. If we have jumped out, we revert to the previous P0 and shrink the step size. At the same time, we also shrink the coefficient of the logarithmic term in the modified objective function. In fact, it will be convenient to give the same value to the step size and the logarithmic coefficient at all times.

Let us see how this works with Problem 11 from Section 6.3 of Helzer's text.

syms x1 x2 x3

c1=126-6*x1-8*x2-7*x3

c2=155-9*x1-8*x2-9*x3

c3=135-10*x1-7*x2-2*x3

c4=132-9*x1-7*x2-5*x3

c5=84-7*x1-2*x2-5*x3

constraints=[c1,c2,c3,c4,c5]

logterm=sum(log(constraints))

object=162*x1+150*x2+128*x3

c1 =

126-6*x1-8*x2-7*x3

c2 =

155-9*x1-8*x2-9*x3

c3 =

135-10*x1-7*x2-2*x3

c4 =

132-9*x1-7*x2-5*x3

c5 =

84-7*x1-2*x2-5*x3

constraints =

[ 126-6*x1-8*x2-7*x3, 155-9*x1-8*x2-9*x3, 135-10*x1-7*x2-2*x3, 132-9*x1-7*x2-5*x3, 84-7*x1-2*x2-5*x3]

logterm =

log(126-6*x1-8*x2-7*x3)+log(155-9*x1-8*x2-9*x3)+log(135-10*x1-7*x2-2*x3)+log(132-9*x1-7*x2-5*x3)+log(84-7*x1-2*x2-5*x3)

object =

162*x1+150*x2+128*x3

We have written a MATLAB function, logbar, which calls a script called barloop and another function called findcrit. logbar takes five arguments which are, in order, the objective function, a row vector whose components are the slack variables (here called constraints), a row vector whose components are the original variables (vars), a starting point (point) and a positive number, which serves as both the step size for findcrit, the critical point finder, and the coefficient of the logarithmic term in the modified objective function. logbar executes the procedure we have outlined above for finding a critical point of the modified objective function until it jumps out of the feasible region. It then returns the previous point. We can then look at the corresponding values of the slack variables. This will tell us when we are close to a vertex of the feasible region. If this is not clear that this is the case, we can reduce the last argument to logbar and evaluate it again. Eventually, it should become apparent that the process is converging to one of the vertices of the feasible region; at that point we can terminate the process and declare the vertex in question to be the solution we are looking for.

vars=[x1,x2,x3]

vars =

[ x1, x2, x3]

We can start at the origin, since it is in the feasible region.

point=[0,0,0]

point =

0 0 0

point=logbar(object, constraints,vars,point,5)

test=subs(constraints,vars,point)

point =

6.2958 8.4314 2.6177

test =

2.4499 7.3269 7.7865 3.2292 9.9777

It is not yet clear that we are near a vertex, so let us try again with a smaller step size and logarithmic term.

point=logbar(object, constraints,vars,point,1)

point =

6.9592 7.1068 3.8837

test=subs(constraints,vars,point)

test =

0.2041 0.5592 7.8929 0.2009 1.6535

Still not clear; once again

point=logbar(object, constraints,vars,point,.1)

point =

6.9920 7.0373 3.9477

test=subs(constraints,vars,point)

test =

0.1156 0.2442 7.9234 0.0723 1.2430

And once more:

point=logbar(object, constraints,vars,point,.02)

test=subs(constraints,vars,point)

point =

7.0003 7.0012 3.9948

test =

0.0253 0.0349 7.9993 0.0152 1.0217

It is now fairly clear that we are approaching the vertex at which the first, second and fourth slack variables will be zero, but let us run one more iteration to be sure.

point=logbar(object, constraints,vars,point,.005)

test=subs(constraints,vars,point)

point =

7.0003 7.0012 3.9948

test =

0.0253 0.0349 7.9993 0.0152 1.0217

Nothing changed; the root finder must have jumped out of the critical region on the first attempt, so we must make the step size even smaller.

point=logbar(object, constraints,vars,point,.001)

test=subs(constraints,vars,point)

point =

7.0000 7.0001 3.9996

test =

0.0015 0.0022 7.9998 0.0009 1.0016

Now it is completely clear.

[x1sol,x2sol,x3sol]=solve(c1,c2,c4)

x1sol =

7

x2sol =

7

x3sol =

4

subs(object,[x1,x2,x3],[x1sol,x2sol,x3sol])

ans =

2696

The maximum of the objective function on the feasible region is 2696, at the point [7,7,4]. Your solution of this problem by the Simplex algorithm should confirm this result.

................
................

#### To fulfill the demand for quickly locating and searching documents.

It is intelligent file search solution for home and business.