@testset "Production Problem" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    c = [4.0, 3.0]
    A1 = [2.0, 1.0, 1.0]
    A2 = [1.0, 2.0, 1.0]
    b1 = 4.0
    b2 = 4.0

    x = MOI.add_variables(optimizer, length(c))

    @test typeof(x[1]) == MOI.VariableIndex

    w, cw = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    cons1 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(A1, [x[1], x[2], y]),
        0.0,
    )
    MOI.add_constraint(optimizer, cons1, MOI.LessThan(b1))

    cons2 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(A2, [x[1], x[2], z]),
        0.0,
    )
    MOI.add_constraint(optimizer, cons2, MOI.LessThan(b2))

    @test cons1.terms[1].coefficient == 2
    @test POI.is_parameter_in_model(optimizer, cons2.terms[3].variable_index)

    obj_func = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([c[1], c[2], 3.0], [x[1], x[2], w]),
        0.0,
    )
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    MOI.get(optimizer, MOI.TerminationStatus())

    MOI.get(optimizer, MOI.PrimalStatus())

    @test isapprox.(
        MOI.get(optimizer, MOI.ObjectiveValue()),
        28 / 3,
        atol = ATOL,
    )
    @test isapprox.(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        4 / 3,
        atol = ATOL,
    )
    @test isapprox.(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        4 / 3,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -5 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -2 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        3.0,
        atol = ATOL,
    )

    MOI.get(optimizer, MOI.VariablePrimal(), w)
    MOI.get(optimizer, MOI.VariablePrimal(), y)
    MOI.get(optimizer, MOI.VariablePrimal(), z)

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(2.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(1.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 2.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) == 1.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) == 1.0

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 13.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [1.0, 1.0]

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -5 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -2 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        3.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(0.0))

    MOI.optimize!(optimizer)

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) == 1.0
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) == 1.0
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 7, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [1.0, 1.0]

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(-5.0))
    MOI.optimize!(optimizer)
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 12.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [3.0, 0.0]

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 0, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -4,
        atol = ATOL,
    )
end

@testset "Production Problem variation on parameters for duals" begin
    optimizer = POI.ParametricOptimizer(GLPK.Optimizer())

    c = [4.0, 3.0]
    A1 = [2.0, 1.0, 3.0]
    A2 = [1.0, 2.0, 0.5]
    b1 = 4.0
    b2 = 4.0

    x = MOI.add_variables(optimizer, length(c))

    @test typeof(x[1]) == MOI.VariableIndex

    w, cw = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    y, cy = MOI.add_constrained_variable(optimizer, POI.Parameter(0))
    z, cz = MOI.add_constrained_variable(optimizer, POI.Parameter(0))

    @test MOI.get(optimizer, MOI.VariablePrimal(), w) == 0

    for x_i in x
        MOI.add_constraint(
            optimizer,
            MOI.SingleVariable(x_i),
            MOI.GreaterThan(0.0),
        )
    end

    cons1 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(A1, [x[1], x[2], y]),
        0.0,
    )
    ci1 = MOI.add_constraint(optimizer, cons1, MOI.LessThan(b1))

    cons2 = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.(A2, [x[1], x[2], z]),
        0.0,
    )
    ci2 = MOI.add_constraint(optimizer, cons2, MOI.LessThan(b2))

    @test cons1.terms[1].coefficient == 2
    @test POI.is_parameter_in_model(optimizer, cons2.terms[3].variable_index)

    obj_func = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([c[1], c[2], 2.0], [x[1], x[2], w]),
        0.0,
    )
    MOI.set(
        optimizer,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        obj_func,
    )
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    MOI.optimize!(optimizer)

    MOI.get(optimizer, MOI.TerminationStatus())

    MOI.get(optimizer, MOI.PrimalStatus())

    @test isapprox(
        MOI.get(optimizer, MOI.ObjectiveValue()),
        28 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[1]),
        4 / 3,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.VariablePrimal(), x[2]),
        4 / 3,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -5,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -2 / 6,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        2.0,
        atol = ATOL,
    )

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        3 * MOI.get(optimizer, MOI.ConstraintDual(), ci1),
        atol = 1e-4,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        0.5 * MOI.get(optimizer, MOI.ConstraintDual(), ci2),
        atol = 1e-4,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(2.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(1.0))
    MOI.set(optimizer, MOI.ConstraintSet(), cz, POI.Parameter(1.0))

    MOI.optimize!(optimizer)

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 7.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [0.0, 1.0]

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -9.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -0.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        2.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cw, POI.Parameter(0.0))

    MOI.optimize!(optimizer)

    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 3.0, atol = ATOL)
    @test MOI.get.(optimizer, MOI.VariablePrimal(), x) == [0.0, 1.0]

    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cy),
        -9.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -0.0,
        atol = ATOL,
    )
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cw),
        2.0,
        atol = ATOL,
    )

    MOI.set(optimizer, MOI.ConstraintSet(), cy, POI.Parameter(-5.0))
    MOI.optimize!(optimizer)
    @test isapprox.(MOI.get(optimizer, MOI.ObjectiveValue()), 14.0, atol = ATOL)
    @test isapprox(
        MOI.get.(optimizer, MOI.VariablePrimal(), x),
        [3.5, 0.0],
        atol = ATOL,
    )

    @test isapprox(MOI.get(optimizer, MOI.ConstraintDual(), cy), 0, atol = ATOL)
    @test isapprox(
        MOI.get(optimizer, MOI.ConstraintDual(), cz),
        -2,
        atol = ATOL,
    )
end
