Programming

Under construction!


On this page you will find programs that I have written, including:

Linear programming:

These versions have been coded as Excel macros so they are easy to use.

The simplex problem is a series of Hilbert matrices. Hilbert matrices are known to be unstable. The code is good for degree 10.

The assignment or transportation problem uses the same simplex code internally:

      Option Explicit
      Sub Simplex()
      Application.ScreenUpdating = False
      Application.Calculation = xlCalculationManual
      Application.EnableEvents = False
      
      Const ZT = 1E-16
      
      Dim Z As Double
      Dim A As Double
      Dim I As Integer
      Dim J As Integer
      Dim M As Integer
      Dim N As Integer
      Dim NM As Integer
      Dim P As Integer
      Dim Q As Integer
      Dim Phase As Integer
      Dim T() As Double
      Dim R() As Double
      Dim B() As Double
      Dim X() As Double

      ' Read model size
      N = ActiveSheet.Cells(1, 1).Value ' Variables (columns)
      M = ActiveSheet.Cells(1, 2).Value ' Constraints (rows)
      NM = N + M                        ' Total columns with slack variables
      ReDim T(1 To M, 1 To NM) As Double
      ReDim R(1 To NM) As Double
      ReDim B(1 To M) As Double
      ReDim S(1 To M) As Double
      ReDim X(1 To N) As Double
        
      ' Read in model
      ' N=Variables, M=Constraints
      ' Maximise:    R(i)
      ' Subject to:  T(j,i) <= B(j)
      '              X(i)   >= 0
        
      For J = 1 To M
        S(J) = 0 ' Used to record basis variables
        For I = 1 To N
          T(J, I) = ActiveSheet.Cells(J + 2, I).Value
        Next I
        B(J) = ActiveSheet.Cells(J + 2, N + 1).Value
      Next J
      For I = 1 To N
        R(I) = ActiveSheet.Cells(2, I).Value
      Next I
      Z = 0 ' System value
    
      ' Add slack variables
      For J = 1 To M
        For I = 1 To M
          T(J, N + I) = 0
        Next I
        R(N + J) = 0
        T(J, N + J) = 1
      Next J
    
      Rem SIMPLEX - Two phase method
      Phase = 1
      While (Phase = 1) Or (Phase = 2)
       
        If (Phase = 1) Then
          Rem Test if feasible - Find smallest B
          Q = 1
          For J = 2 To M
            If (B(Q) > B(J)) Then Q = J
          Next J
          If (B(Q) >= 0) Then
            Rem Feasible but not necessarily optimal
            Phase = Phase + 1
          Else
            Rem Not feasible - Select special pivot - Find smallest ratio
            P = 0
            For I = 1 To N + M
              If (T(Q, I) < 0) Then
                If (P = 0) Then P = I
                If R(P) / T(Q, P) > (R(I) / T(Q, I)) Then P = I
              End If
            Next I
            If (P = 0) Then
              Rem Infeasible
              Phase = -2
            End If
          End If
        End If
        
        If Phase = 2 Then
          Rem Test if optimal - Find largest R
          P = 1
          For I = 2 To N + M
            If (R(P) < R(I)) Then P = I
          Next I
          If (R(P) <= 0) Then
            Rem Optimal
            Phase = Phase + 1
          Else
            Rem Not optimal - Select pivot - Find smallest ratio
            Q = 0
            For J = 1 To M
              If (T(J, P) > 0) Then
                If (Q = 0) Then Q = J
                If (B(Q) / T(Q, P) > B(J) / T(J, P)) Then Q = J
              End If
            Next J
            If (Q = 0) Then
              Rem Unbounded
              Phase = -1
            End If
          End If
        End If
        
        Rem Row operations - Test for near zero here
        If (Phase = 1) Or (Phase = 2) Then
          S(Q) = P
          A = T(Q, P)
          For I = 1 To N + M
            T(Q, I) = T(Q, I) / A
            If (Abs(T(Q, I)) <= ZT) Then T(Q, I) = 0
          Next I
          T(Q, P) = 1
          B(Q) = B(Q) / A
          If (Abs(B(Q)) <= ZT) Then B(Q) = 0
          For J = 1 To M
            If (J <> Q) Then
              For I = 1 To N + M
               If (I <> P) Then
                  T(J, I) = T(J, I) - T(J, P) * T(Q, I)
                  If (Abs(T(J, I)) <= ZT) Then T(J, I) = 0
                End If
              Next I
              B(J) = B(J) - T(J, P) * B(Q)
              If (Abs(B(J)) <= ZT) Then B(J) = 0
            End If
          Next J
          For J = 1 To M
            T(J, P) = 0
          Next J
          T(Q, P) = 1
          A = R(P)
          For I = 1 To N + M
            R(I) = R(I) - A * T(Q, I)
            If (Abs(R(I)) <= ZT) Then R(I) = 0
          Next I
          R(P) = 0
          Z = Z + A * B(Q)
        End If
    
      Wend
  
      If (Phase = 3) Then
        ' Decode solution
        For I = 1 To N
          X(I) = 0
          For J = 1 To M
            If (S(J) = I) Then X(I) = B(J)
          Next J
          ActiveSheet.Cells(3 + M, I).Value = X(I)
        Next I
        ActiveSheet.Cells(3 + M, N + 1).Value = Z
      Else
        For I = 1 To N
          ActiveSheet.Cells(3 + M, I).Value = ""
        Next I
        ActiveSheet.Cells(3 + M, N + 1).Value = "No Solution"
      End If

      Application.ScreenUpdating = True
      Application.Calculation = xlCalculationAutomatic
      Application.EnableEvents = True
      End Sub
    

You can contact me at Master@AlanCooper.id.au if you feel the need.


Home