Generalized linear programming solver

LpSolver.kt

import kotlin.math.abs

/*
 * Generalized linear programming solver.
 *
 * What changed relative to the original single-purpose simplex:
 *
 *   1. Constraints may be <=, >=, or ==, with rhs of any sign.
 *      A two-phase simplex (phase 1 minimizes artificial variables)
 *      finds a feasible starting basis instead of assuming slack-only.
 *   2. The basis is tracked explicitly in an IntArray. Solution
 *      extraction reads the basis directly; no unit-column detection.
 *   3. Dantzig pivoting with a Bland's-rule fallback after an
 *      iteration budget, so degenerate problems cannot cycle.
 *   4. Outcomes form a closed ADT: Optimal / Infeasible / Unbounded.
 *   5. A modeling layer (LpModel) handles variable lower/upper bounds
 *      by shifting internally, generalizing the manual
 *      "actual = min + extra" transformation.
 */

// ---------------------------------------------------------------------------
// Result ADT
// ---------------------------------------------------------------------------

sealed interface LpResult {
    data class Optimal(
        val values: Map<LpVariable, Double>,
        val objective: Double
    ) : LpResult

    data object Infeasible : LpResult
    data object Unbounded : LpResult
}

enum class Relation { LE, GE, EQ }

// ---------------------------------------------------------------------------
// Modeling layer
// ---------------------------------------------------------------------------

class LpVariable internal constructor(
    val name: String,
    val lower: Double,
    val upper: Double
) {
    override fun toString() = name
}

class LinExpr internal constructor(
    internal val coeffs: Map<LpVariable, Double>
) {
    operator fun plus(other: LinExpr): LinExpr {
        val merged = coeffs.toMutableMap()
        for ((v, c) in other.coeffs) merged.merge(v, c, Double::plus)
        return LinExpr(merged)
    }

    operator fun plus(v: LpVariable) = this + v.expr()
    operator fun minus(other: LinExpr) = this + (-1.0 * other)
    operator fun minus(v: LpVariable) = this + (-1.0 * v)

    infix fun le(rhs: Double) = LpConstraint(this, Relation.LE, rhs)
    infix fun ge(rhs: Double) = LpConstraint(this, Relation.GE, rhs)
    infix fun eq(rhs: Double) = LpConstraint(this, Relation.EQ, rhs)

    fun evaluate(values: Map<LpVariable, Double>): Double =
        coeffs.entries.sumOf { (v, c) -> c * values.getValue(v) }
}

fun LpVariable.expr() = LinExpr(mapOf(this to 1.0))

operator fun Double.times(v: LpVariable) = LinExpr(mapOf(v to this))
operator fun Double.times(e: LinExpr) =
    LinExpr(e.coeffs.mapValues { (_, c) -> this * c })

operator fun LpVariable.plus(other: LpVariable) = expr() + other
operator fun LpVariable.plus(e: LinExpr) = expr() + e
operator fun LpVariable.minus(other: LpVariable) = expr() - other
operator fun LpVariable.minus(e: LinExpr) = expr() - e

infix fun LpVariable.le(rhs: Double) = expr() le rhs
infix fun LpVariable.ge(rhs: Double) = expr() ge rhs
infix fun LpVariable.eq(rhs: Double) = expr() eq rhs

class LpConstraint internal constructor(
    internal val expr: LinExpr,
    internal val relation: Relation,
    internal val rhs: Double
)

class LpModel {
    private val variables = mutableListOf<LpVariable>()
    private val constraints = mutableListOf<LpConstraint>()
    private var objective: LinExpr? = null
    private var maximize = true

    fun variable(
        name: String,
        lower: Double = 0.0,
        upper: Double = Double.POSITIVE_INFINITY
    ): LpVariable {
        require(lower.isFinite()) {
            "$name: free variables are not supported; give a finite lower bound " +
                "(or split into positive and negative parts)"
        }
        require(lower <= upper) { "$name: lower bound exceeds upper bound" }
        return LpVariable(name, lower, upper).also { variables += it }
    }

    fun maximize(expr: LinExpr) { objective = expr; maximize = true }
    fun minimize(expr: LinExpr) { objective = expr; maximize = false }

    fun add(constraint: LpConstraint) { constraints += constraint }

    fun solve(): LpResult {
        val obj = requireNotNull(objective) { "No objective set" }
        val index = variables.withIndex().associate { (i, v) -> v to i }
        val n = variables.size

        // Internal space: x = lower + x', with x' >= 0.
        fun rowOf(expr: LinExpr): DoubleArray {
            val row = DoubleArray(n)
            for ((v, c) in expr.coeffs) row[index.getValue(v)] += c
            return row
        }

        fun shiftOf(expr: LinExpr): Double =
            expr.coeffs.entries.sumOf { (v, c) -> c * v.lower }

        val rows = mutableListOf<DoubleArray>()
        val relations = mutableListOf<Relation>()
        val rhs = mutableListOf<Double>()

        for (c in constraints) {
            rows += rowOf(c.expr)
            relations += c.relation
            rhs += c.rhs - shiftOf(c.expr)
        }

        // Finite upper bounds become x' <= upper - lower.
        for ((i, v) in variables.withIndex()) {
            if (v.upper.isFinite()) {
                rows += DoubleArray(n).also { it[i] = 1.0 }
                relations += Relation.LE
                rhs += v.upper - v.lower
            }
        }

        val solver = TwoPhaseSimplex(
            maximize = maximize,
            objective = rowOf(obj),
            rows = rows.toTypedArray(),
            relations = relations.toTypedArray(),
            rhs = rhs.toDoubleArray()
        )

        return when (val outcome = solver.solve()) {
            is TwoPhaseSimplex.Outcome.Infeasible -> LpResult.Infeasible
            is TwoPhaseSimplex.Outcome.Unbounded -> LpResult.Unbounded
            is TwoPhaseSimplex.Outcome.Optimal -> {
                val values = variables.withIndex().associate { (i, v) ->
                    v to v.lower + outcome.x[i]
                }
                LpResult.Optimal(values, obj.evaluate(values))
            }
        }
    }
}

// ---------------------------------------------------------------------------
// Solver core
// ---------------------------------------------------------------------------

class TwoPhaseSimplex(
    private val maximize: Boolean,
    private val objective: DoubleArray,
    rows: Array<DoubleArray>,
    relations: Array<Relation>,
    rhs: DoubleArray
) {
    sealed interface Outcome {
        data class Optimal(val x: DoubleArray) : Outcome
        data object Infeasible : Outcome
        data object Unbounded : Outcome
    }

    private val m = rows.size
    private val n = objective.size

    // Normalized copies: rhs >= 0 (negating a row flips its relation).
    private val rows: Array<DoubleArray>
    private val relations: Array<Relation>
    private val rhs: DoubleArray

    init {
        require(rows.all { it.size == n })
        require(relations.size == m && rhs.size == m)

        this.rows = Array(m) { rows[it].copyOf() }
        this.relations = relations.copyOf()
        this.rhs = rhs.copyOf()

        for (i in 0 until m) {
            if (this.rhs[i] < 0) {
                for (j in 0 until n) this.rows[i][j] = -this.rows[i][j]
                this.rhs[i] = -this.rhs[i]
                this.relations[i] = when (this.relations[i]) {
                    Relation.LE -> Relation.GE
                    Relation.GE -> Relation.LE
                    Relation.EQ -> Relation.EQ
                }
            }
        }
    }

    // Column layout: [structural | slack and surplus | artificial | rhs]
    private val slackCount = relations.count { it != Relation.EQ }
    private val artificialCount = relations.count { it != Relation.LE }
    private val cols = n + slackCount + artificialCount
    private val tableau = Array(m + 1) { DoubleArray(cols + 1) }
    private val basis = IntArray(m)
    private val artificialColumns = mutableSetOf<Int>()

    private companion object {
        const val EPS = 1e-9
        const val FEASIBILITY_TOL = 1e-7
    }

    fun solve(): Outcome {
        buildTableau()

        if (artificialColumns.isNotEmpty()) {
            when (phase1()) {
                Phase1Result.INFEASIBLE -> return Outcome.Infeasible
                Phase1Result.FEASIBLE -> Unit
            }
        }

        return when (phase2()) {
            SimplexStatus.UNBOUNDED -> Outcome.Unbounded
            SimplexStatus.OPTIMAL -> Outcome.Optimal(extract())
        }
    }

    private fun buildTableau() {
        var slack = n
        var artificial = n + slackCount

        for (i in 0 until m) {
            for (j in 0 until n) tableau[i][j] = rows[i][j]
            tableau[i][cols] = rhs[i]

            when (relations[i]) {
                Relation.LE -> {
                    tableau[i][slack] = 1.0
                    basis[i] = slack
                    slack++
                }
                Relation.GE -> {
                    tableau[i][slack] = -1.0
                    slack++
                    tableau[i][artificial] = 1.0
                    basis[i] = artificial
                    artificialColumns += artificial
                    artificial++
                }
                Relation.EQ -> {
                    tableau[i][artificial] = 1.0
                    basis[i] = artificial
                    artificialColumns += artificial
                    artificial++
                }
            }
        }
    }

    private enum class Phase1Result { FEASIBLE, INFEASIBLE }
    private enum class SimplexStatus { OPTIMAL, UNBOUNDED }

    private fun phase1(): Phase1Result {
        // Minimize the sum of artificials, i.e. maximize its negation.
        // With the negated-coefficient convention the objective row
        // holds +1 for each artificial column.
        tableau[m].fill(0.0)
        for (a in artificialColumns) tableau[m][a] = 1.0
        priceOutBasicColumns()

        when (runSimplex(allowArtificials = true)) {
            SimplexStatus.UNBOUNDED ->
                error("Phase 1 objective is bounded by construction")
            SimplexStatus.OPTIMAL -> Unit
        }

        val artificialSum = -tableau[m][cols]
        if (artificialSum > FEASIBILITY_TOL) return Phase1Result.INFEASIBLE

        driveOutArtificials()
        return Phase1Result.FEASIBLE
    }

    private fun phase2(): SimplexStatus {
        val sign = if (maximize) 1.0 else -1.0
        tableau[m].fill(0.0)
        for (j in 0 until n) tableau[m][j] = -sign * objective[j]
        priceOutBasicColumns()

        return runSimplex(allowArtificials = false)
    }

    /** The objective row must have zeros in basic columns. */
    private fun priceOutBasicColumns() {
        for (i in 0 until m) {
            val factor = tableau[m][basis[i]]
            if (abs(factor) > EPS) {
                for (j in 0..cols) tableau[m][j] -= factor * tableau[i][j]
            }
        }
    }

    /**
     * Artificials still basic at level zero after phase 1 are pivoted out
     * on any usable structural or slack column; a row with no such column
     * is redundant and stays parked on its artificial at value zero.
     */
    private fun driveOutArtificials() {
        for (i in 0 until m) {
            if (basis[i] !in artificialColumns) continue

            val col = (0 until cols).firstOrNull { j ->
                j !in artificialColumns && abs(tableau[i][j]) > EPS
            }

            if (col != null) {
                pivot(i, col)
                basis[i] = col
            }
        }
    }

    private fun runSimplex(allowArtificials: Boolean): SimplexStatus {
        // Dantzig's rule is fast but can cycle on degenerate problems;
        // past the iteration budget we switch to Bland's rule, which
        // is guaranteed to terminate.
        val dantzigBudget = 50 * (m + cols)
        var iterations = 0

        while (true) {
            iterations++
            val bland = iterations > dantzigBudget

            val pivotCol = enteringColumn(allowArtificials, bland) ?: return SimplexStatus.OPTIMAL
            val pivotRow = leavingRow(pivotCol) ?: return SimplexStatus.UNBOUNDED

            pivot(pivotRow, pivotCol)
            basis[pivotRow] = pivotCol
        }
    }

    private fun enteringColumn(allowArtificials: Boolean, bland: Boolean): Int? {
        var col = -1
        var mostNegative = -EPS

        for (j in 0 until cols) {
            if (!allowArtificials && j in artificialColumns) continue

            val value = tableau[m][j]
            if (value < mostNegative) {
                if (bland) return j // first eligible column
                mostNegative = value
                col = j
            }
        }

        return col.takeIf { it != -1 }
    }

    private fun leavingRow(pivotCol: Int): Int? {
        var bestRow: Int? = null
        var bestRatio = Double.POSITIVE_INFINITY

        for (i in 0 until m) {
            val coefficient = tableau[i][pivotCol]
            if (coefficient <= EPS) continue

            val ratio = tableau[i][cols] / coefficient
            val strictlyBetter = ratio < bestRatio - EPS
            val tieOnLowerBasis = abs(ratio - bestRatio) <= EPS &&
                bestRow != null && basis[i] < basis[bestRow]

            if (strictlyBetter || tieOnLowerBasis || bestRow == null && ratio < bestRatio) {
                bestRatio = ratio
                bestRow = i
            }
        }

        return bestRow
    }

    private fun pivot(pivotRow: Int, pivotCol: Int) {
        val pivotValue = tableau[pivotRow][pivotCol]
        for (j in 0..cols) tableau[pivotRow][j] /= pivotValue

        for (i in 0..m) {
            if (i == pivotRow) continue
            val factor = tableau[i][pivotCol]
            if (abs(factor) > 0.0) {
                for (j in 0..cols) tableau[i][j] -= factor * tableau[pivotRow][j]
            }
        }
    }

    private fun extract(): DoubleArray {
        val x = DoubleArray(n)
        for (i in 0 until m) {
            if (basis[i] < n) x[basis[i]] = tableau[i][cols]
        }
        return x
    }
}

usage.kt

fun main() {
    val model = LpModel()

    val x = model.variable("x", lower = 0.0, upper = 10.0)
    val y = model.variable("y", lower = 0.0)  // upper defaults to +infinity

    model.add(x + y le 12.0)
    model.add(2.0 * x - y ge 0.0)

    model.maximize(3.0 * x + 5.0 * y)

    when (val result = model.solve()) {
        is LpResult.Optimal -> {
            println("x = ${result.values.getValue(x)}")
            println("y = ${result.values.getValue(y)}")
            println("objective = ${result.objective}")
        }
        LpResult.Infeasible -> println("no solution satisfies the constraints")
        LpResult.Unbounded -> println("objective grows without limit")
    }
}