 ```! LINGO MODEL for interior point/barrier primal-dual method. (IntPtAlg.lg4) Original version by Kipp Martin. See article by Marsten, et. al.(1990) in INTERFACES, vol. 20, no. 4.; ! Keywords: Interior Point Method, Barrier Method; SETS: ROW : B, Y, DY; COL : C, DX, X, DZ, Z; MATRIX( ROW, COL) : A; ENDSETS ! Solve the optimization problem Min CX s.t. AX = B, X >= 0, by solving approximately the system: 1) AX = B, (Primal conditions) 2) A'Y - C = Z,(Dual conditions) 3) X >= 0, (Non-negativity) 4) Z >= 0 5) ZX = 0. (Complementarity) We start with a solution to the linear constraints, (1-4), and then use iterative approximation to satisfy (5); DATA: COL = A__ C__ D__ S_1 S_2 S_3; B = 60, 50, 120; ! RHS; C = -20, -30, -47, 0, 0, 0; ! Objective coefficients, cost minimize form; A = 1, 0, 1, 1, 0, 0, ! Coefficient matrix for the; 0, 1, 1, 0, 1, 0, ! Astro/Cosmo; 1, 2, 3, 0, 0, 1; ! Problem in Schrage; X = 1 1 1 58 48 114; ! Initial feasible primal soln; Z = 40 60 103 30 30 30; ! Initial feasible reduced costs; Y = 30 30 30; ! Initial dual prices; TOLZR = .0000001; ! Zero tolerance; TOLGAP = .00000005; ! Termination tolerance; ENDDATA SUBMODEL DIRECTION: ! Given a current point: X : feasible (nonnegative) primal point, Y : feasible dual prices, Z : feasible (nonnegative) reduced costs, Choose a direction: DX, DY and DZ that tries to better satisfy the complementary slackness constraint; ! The Primal equations: A*( X + DX) = B; @FOR( ROW( I): [P] @SUM( COL( J) : A( I, J)*( X( J) + DX( J))) = B( I); ); ! The Dual equations: A*( Y + DY) -(Z + DZ) = -C; @FOR( COL( J): [D] @SUM( ROW( I): A( I, J) * ( Y( I) + DY( I))) - (Z( J) + DZ( J)) = -C( J) ); ! The Complementary slackness equation: ( Z + DZ) * ( X + DX) = 0, or Z*X + Z*DX + DZ*X + DZ*DX = 0, or the linear approximation Z*X + Z*DX + DZ*X ~ 0, ; ! and the linear approximation to it; @FOR( COL( J): [CM] ( Z( J) + DZ( J)) * X( J) + Z( J) * DX( J) = MU ); ! Allow for movements in the negative direction; @FOR( ROW: @FREE( DY)); @FOR( COL: @FREE( DX)); @FOR( COL: @FREE( DZ)); ENDSUBMODEL CALC: N = @SIZE( COL); ! Turn off all messages; @SET('TERSEO',3); ! Set to <= 2 for debugging; maxiter = 30; ! Seldom need more than this # iterations; more = 1; ! Set to 0 when done; iter = 0; @WHILE( iter #LE# maxiter #AND# more: iter = iter + 1; gap = @SUM( COL: X * Z); ! Duality gap; @WRITE(iter,') Current gap= ', gap,@NEWLINE(1)); ! Marsten's rule for error target, MU; MU = gap / ( N * N); @SOLVE( DIRECTION); ! Get the move direction; ! Calculate step size in the move direction so that X and Z remain > 0; ! Do not let any variable get too close to zero; ALPHAP = .995 * ( @MIN(COL( J) | DX( J) #LT# -TOLZR: -X(J) / DX(J))); ALPHAD = .995 * ( @MIN(COL( J) | DZ( J) #LT# -TOLZR: -Z(J) / DZ(J))); ! Move to the new point; @FOR( COL( J): X( J) = X( J) + ALPHAP * DX( J)); @FOR( COL( J): Z( J) = Z( J) + ALPHAD * DZ( J)); @FOR( ROW( I): Y( I) = Y( I) + ALPHAD * DY( I)); ! Calculate primal and dual objectives; POBJ = -@SUM( COL: X * C); DOBJ = @SUM( ROW: Y * B); gap = DOBJ - POBJ ; ! Write out details of new point; @WRITE(' New Primal, Dual, Gap= ', @FORMAT(POBJ,'12.4f'),' ',@FORMAT(DOBJ,'12.4f'),' ',gap,@NEWLINE(1)); @FOR( COL(j): @WRITE(' ',COL(j))); @WRITE(@NEWLINE(1),'X= '); @FOR( COL(j): @WRITE( @FORMAT(X(j),'10.3f'),' '); ); @WRITE(@NEWLINE(1)); @WRITE('Z= '); @FOR( COL(j): @WRITE( @FORMAT(Z(j),'10.3f'),' '); ); @WRITE(@NEWLINE(1)); @WRITE('Y= '); @FOR( ROW(i): @WRITE( @FORMAT(Y(i),'10.3f'),' '); ); @WRITE(@NEWLINE(1)); ! Exit if gap is small enough; @IFC( gap #LE# TOLGAP *( 1 + @ABS( @SUM( COL: C * X))): more = 0); ); ENDCALC ```