How do I identify incorrect coding to get the correct result from my R linear programming optimisation model?

Question

I have a linear programming task that I am struggling with. I have formulated the LP problem and have generated some code that appears to function; however, I am suspicious of the outputs and do not believe it to be the correct answer to the question.

**The problem is this: ** A factory makes three products called Spring, Autumn, and Winter, from three materials containing Cotton, Wool and Silk.

            Sale Price    Production Cost
    Spring  60            5
    Autumn  55            3
    Winter  60            5

            Purchase Price 
    Cotton  30
    Wool    45
    Silk    50
    
            Demand  Min Cotton Proportion   Min Wool Proportion
    Spring  3300    .55                    .30
    Autumn  3600    .45                    .4
    Winter  4000    .3                     .5

My Formulation of the problem is this:

Decision variables:
xij ≥ 0 denotes # of tons of
product_j ϵ {1 = Spring, 2 = Autumn, 3 = Winter}
produced from
material_i ϵ {c = Cotton, w = Wool, s = Silk}

So let:
1c, 1w, 1s, = cotton, wool and silk demand for spring
2c, 2w, 2s = cotton, wool and silk demand for autumn
3c, 3w, 2s = cotton, wool and silk demand for winter

Objective Function:
Maximize profit
Profit = Sales Price - Production Cost - Purchase Price

MAX ((55 * 1c) + (55 * 1w) + (55 * 1s) + (52 * 2c) + (52 * 2w) + (52 * 2s) + (55 * 3c) + (55 * 3w) + (55 * 3W)) - ((30 * 1c) - (30 * 2c) - (30 * 3c) - (45 * 1w) - (45 * 2w) - (45 * 3w) - (50 * 1s) - (50 * 2s) - (50 * 3s))

Constraints:
Spring Cotton Proportion + Spring Wool Proportion + Spring silk Proportion <= 3300 or 1c + 1w + 1s <= 3300
Autumn Cotton Proportion + Autumn Wool Proportion + Autumn silk Proportion <= 3600 or 2c + 2w + 2s <= 3600
Winter Cotton Proportion + Winter Wool Proportion + Winter silk Proportion <= 400 or 3c + 3w + 3s <= 4000

Spring Cotton Proportion >= 55% or 1c >= .55
Autumn Cotton Proportion >= 45% or 2c >= .45
Winter Cotton Proportion >= 30% or 3c >= .30
Spring Wool Proportion >= 30% or 1w >= .30
Autumn Wool Proportion >= 40% or 2w >= .40
Winter Wool Proportion >= 50% or 3w >= .50

All variables >= 0

# Define the LP model
model2 <- make.lp(0, 9)  # 3 rows for demand constraints, 9 columns for decision variables

# Set the objective function coefficients
obj.coefs <- c(55, 55, 55, 52, 52, 52, 55, 55, 55) - c(30, 30, 30, 45, 45, 45, 50, 50, 50)
lp.control(model2, sense = "max")  # Set the optimization sense to maximize
$anti.degen
[1] "fixedvars" "stalling" 

$basis.crash
[1] "none"

$bb.depthlimit
[1] -50

$bb.floorfirst
[1] "automatic"

$bb.rule
[1] "pseudononint" "greedy"       "dynamic"      "rcostfixing" 

$break.at.first
[1] FALSE

$break.at.value
[1] 1e+30

$epsilon
      epsb       epsd      epsel     epsint epsperturb   epspivot 
     1e-10      1e-09      1e-12      1e-07      1e-05      2e-07 

$improve
[1] "dualfeas" "thetagap"

$infinite
[1] 1e+30

$maxpivot
[1] 250

$mip.gap
absolute relative 
   1e-11    1e-11 

$negrange
[1] -1e+06

$obj.in.basis
[1] TRUE

$pivoting
[1] "devex"    "adaptive"

$presolve
[1] "none"

$scalelimit
[1] 5

$scaling
[1] "geometric"   "equilibrate" "integers"   

$sense
[1] "maximize"

$simplextype
[1] "dual"   "primal"

$timeout
[1] 0

$verbose
[1] "neutral"

set.objfn(model2, obj.coefs)

# Add the constraints
add.constraint(model2, c(1, 1, 1, 0, 0, 0, 0, 0, 0), "<=", 3300)  # Demand constraint for Spring (Cotton, Wool, Silk)
add.constraint(model2, c(0, 0, 0, 1,1,1, 0, 0, 0), "<=", 3600)  # Demand constraint for Autumn (Cotton, Wool, Silk)
add.constraint(model2, c(0, 0, 0, 0, 0, 0, 1,1,1), "<=", 4000)  # Demand constraint for Winter (Cotton, Wool, Silk)
# Set the cotton proportion constraints
add.constraint(model2, c(0.55, 0, 0, 0, 0, 0, 0, 0, 0), ">=", 0)  # Spring cotton proportion
add.constraint(model2, c(0, 0, 0, 0.45, 0, 0, 0, 0, 0), ">=", 0)  # Autumn cotton proportion
add.constraint(model2, c(0, 0, 0, 0, 0, 0, 0.30, 0, 0), ">=", 0)  # Winter cotton proportion

# Set the wool proportion constraints
add.constraint(model2, c(0, 0.30, 0, 0, 0, 0, 0, 0, 0), ">=", 0)  # Spring wool proportion
add.constraint(model2, c(0, 0, 0, 0, 0.40, 0, 0, 0, 0), ">=", 0)  # Autumn wool proportion
add.constraint(model2, c(0, 0, 0, 0, 0, 0, 0, 0.50, 0), ">=", 0)  # Winter wool proportion
add.constraint(model2, rep(0, 9), ">", 0)  # Non-negativity constraint for all variables

# Solve the LP model
solve(model2)
[1] 0
 
# Get the optimal profit
optimal_profit <- get.objective(model2)

# Get the optimal values of the decision variables
optimal_values <- get.variables(model2)

# Display the results
cat("Optimal Profit: $", optimal_profit, "\n")
Optimal Profit: $ 127700 
cat("Optimal Values of Decision Variables:\n")
Optimal Values of Decision Variables:
cat("Spring (Cotton):", optimal_values[1], "tons\n")
Spring (Cotton): 3300 tons
cat("Spring (Wool):", optimal_values[2], "tons\n")
Spring (Wool): 0 tons
cat("Spring (Silk):", optimal_values[3], "tons\n")
Spring (Silk): 0 tons
cat("Autumn (Cotton):", optimal_values[4], "tons\n")
Autumn (Cotton): 3600 tons
cat("Autumn (Wool):", optimal_values[5], "tons\n")
Autumn (Wool): 0 tons
cat("Autumn (Silk):", optimal_values[6], "tons\n")
Autumn (Silk): 0 tons
cat("Winter (Cotton):", optimal_values[7], "tons\n")
Winter (Cotton): 4000 tons
cat("Winter (Wool):", optimal_values[8], "tons\n")
Winter (Wool): 0 tons
cat("Winter (Silk):", optimal_values[9], "tons\n")
Winter (Silk): 0 tons

Answer 1

Your code does not implement your stated constraints around material proportion minima, which is why your output looks incorrect. Once you do, your output should look like

import pandas as pd
import pulp

products = pd.DataFrame(
    index=pd.Index(name='product', data=('Spring', 'Autumn', 'Winter')),
    data={
        'sale_price': (60, 55, 60),
        'production_cost': (5, 3, 5),
        'demand': (3300, 3600, 4000),
    }
)
materials = pd.Series(
    name='purchase_price',
    index=pd.Index(name='material', data=('Cotton', 'Wool', 'Silk')),
    data=(30, 45, 50),
)
proportions = pd.DataFrame(
    index=products.index,
    columns=pd.Index(name='material', data=('Cotton', 'Wool')),
    data=(
        (.55, .30),
        (.45, .40),
        (.30, .50),
    ),
)

print(products, '\n')
print(proportions, '\n')
print(materials, '\n')


def make_vars(row: pd.Series) -> pulp.LpVariable:
    return pulp.LpVariable(
        name=f'{row["product"]}_{row.material}',
        cat=pulp.LpContinuous, lowBound=0)


combinations = pd.MultiIndex.from_product((products.index, materials.index))
production = combinations.to_frame().apply(make_vars, axis=1)
products['production'] = production.unstack(level='material').apply(pulp.lpSum, axis=1)

prob = pulp.LpProblem(name='clothing_production', sense=pulp.LpMaximize)
profit = pulp.LpAffineExpression()

for product, row in products.iterrows():
    prob.addConstraint(
        name=f'{product}_demand',
        constraint=row.production <= row.demand,
    )
    profit += (row.sale_price - row.production_cost)*row.production

for product, row in proportions.iterrows():
    for material, proportion in row.items():
        prob.addConstraint(
            name=f'{product}_{material}_mix',
            constraint=production[(product, material)] >= proportion * products.loc[product, 'production']
        )

for material, outputs in production.groupby(level='material'):
    profit -= pulp.lpSum(outputs) * materials[material]

prob.objective = profit
print(prob)
prob.solve()
assert prob.status == pulp.LpStatusOptimal

production = production.unstack(level='material').applymap(pulp.LpVariable.value)
print('Production')
print(production)
print()

print('Material ratio')
print(production.div(production.sum(axis=1), axis=0))
print()

         sale_price  production_cost  demand
product                                     
Spring           60                5    3300
Autumn           55                3    3600
Winter           60                5    4000 

material  Cotton  Wool
product               
Spring      0.55   0.3
Autumn      0.45   0.4
Winter      0.30   0.5 

material
Cotton    30
Wool      45
Silk      50
Name: purchase_price, dtype: int64 

clothing_production:
MAXIMIZE
22*Autumn_Cotton + 2*Autumn_Silk + 7*Autumn_Wool + 25*Spring_Cotton + 5*Spring_Silk + 10*Spring_Wool + 25*Winter_Cotton + 5*Winter_Silk + 10*Winter_Wool + 0
SUBJECT TO
Spring_demand: Spring_Cotton + Spring_Silk + Spring_Wool <= 3300

Autumn_demand: Autumn_Cotton + Autumn_Silk + Autumn_Wool <= 3600

Winter_demand: Winter_Cotton + Winter_Silk + Winter_Wool <= 4000

Spring_Cotton_mix: 0.45 Spring_Cotton - 0.55 Spring_Silk - 0.55 Spring_Wool
 >= 0

Spring_Wool_mix: - 0.3 Spring_Cotton - 0.3 Spring_Silk + 0.7 Spring_Wool >= 0

Autumn_Cotton_mix: 0.55 Autumn_Cotton - 0.45 Autumn_Silk - 0.45 Autumn_Wool
 >= 0

Autumn_Wool_mix: - 0.4 Autumn_Cotton - 0.4 Autumn_Silk + 0.6 Autumn_Wool >= 0

Winter_Cotton_mix: 0.7 Winter_Cotton - 0.3 Winter_Silk - 0.3 Winter_Wool >= 0

Winter_Wool_mix: - 0.5 Winter_Cotton - 0.5 Winter_Silk + 0.5 Winter_Wool >= 0

VARIABLES
Autumn_Cotton Continuous
Autumn_Silk Continuous
Autumn_Wool Continuous
Spring_Cotton Continuous
Spring_Silk Continuous
Spring_Wool Continuous
Winter_Cotton Continuous
Winter_Silk Continuous
Winter_Wool Continuous

Optimal - objective value 195250
Optimal objective 195250 - 6 iterations time 0.002
Option for printingOptions changed from normal to all
Total time (CPU seconds):       0.00   (Wallclock seconds):       0.00

Production
material  Cotton  Silk    Wool
product                       
Autumn    2160.0   0.0  1440.0
Spring    2310.0   0.0   990.0
Winter    2000.0   0.0  2000.0

Material ratio
material  Cotton  Silk  Wool
product                     
Autumn       0.6   0.0   0.4
Spring       0.7   0.0   0.3
Winter       0.5   0.0   0.5

I leave conversion to R as an exercise to the reader.