Self_Assess_Solns.ipynb

Kernel: Python 3 (system-wide)

Lesson 02 - Self-Assessment Solutions

Self Assessment: Investment Allocation

from pyomo.environ import *

model = ConcreteModel()
model.friend1 = Var(domain=NonNegativeReals)
model.friend2 = Var(domain=NonNegativeReals)

# Objective
model.profit = Objective(expr = 9000 * model.friend1 + 9000 * model.friend2,sense=maximize)

# Constraints
model.fraction1 = Constraint( expr =     1 * model.friend1 +    0 * model.friend2 <= 1)
model.fraction2 = Constraint( expr =     0 * model.friend1 +    1 * model.friend2 <= 1)
model.money =     Constraint( expr = 10000 * model.friend1 + 8000 * model.friend2 <= 12000)
model.workhours = Constraint( expr =   400 * model.friend1 +  500 * model.friend2 <= 600)

# Solve
solver = SolverFactory('glpk')
solver.solve(model)

model.profit.pprint()
model.fraction1.pprint()
model.fraction2.pprint()
model.money.pprint()
model.workhours.pprint()

# display solution
print(f"Maximum Profit = ${model.profit():,.2f}")
print(f"One = {model.friend1():.3f}")
print(f"Two= {model.friend2():.3f}")

profit : Size=1, Index=None, Active=True
    Key  : Active : Sense    : Expression
    None :   True : maximize : 9000*friend1 + 9000*friend2
fraction1 : Size=1, Index=None, Active=True
    Key  : Lower : Body                : Upper : Active
    None :  -Inf : friend1 + 0*friend2 :   1.0 :   True
fraction2 : Size=1, Index=None, Active=True
    Key  : Lower : Body                : Upper : Active
    None :  -Inf : 0*friend1 + friend2 :   1.0 :   True
money : Size=1, Index=None, Active=True
    Key  : Lower : Body                         : Upper   : Active
    None :  -Inf : 10000*friend1 + 8000*friend2 : 12000.0 :   True
workhours : Size=1, Index=None, Active=True
    Key  : Lower : Body                      : Upper : Active
    None :  -Inf : 400*friend1 + 500*friend2 : 600.0 :   True
Maximum Profit = $12,000.00
One = 0.667
Two= 0.667

from pyomo.environ import *

### PROBLEM DATA ###

# Load Data
friends = ['Friend1', 'Friend2']
constraints = ['Fraction1','Fraction2','Money','Work_Hours']
profits = [9000, 9000]
rhs = [1,1,12000,600] # right hand side
coef = [[1,0], 
        [0,1], 
        [10000,8000], 
        [400,500]] # times [f1, f2]

# parse into dictionaries
profit_rate = dict( zip( friends, profits) )
constraint_rhs = dict( zip( constraints, rhs) )
constraint_coef = { c: {f: coef[i][j] for j,f in enumerate(friends)} for i,c in enumerate(constraints)}

### MODEL CONSTRUCTION ###

# Declaration
model = ConcreteModel()

# Decision Variables
model.invest_frac = Var(friends, domain=NonNegativeReals)

# Objective 9000 * friend1 + 9000 * friend2
model.profit = Objective(expr=sum(profit_rate[f] * model.invest_frac[f]
                               for f in friends),
                      sense=maximize)

# model.fraction1 = Constraint( expr =     1 * model.friend1 +    0 * model.friend2 <= 1)
# model.fraction2 = Constraint( expr =     0 * model.friend1 +    1 * model.friend2 <= 1)
# model.money =     Constraint( expr = 10000 * model.friend1 + 8000 * model.friend2 <= 12000)
# model.workhours = Constraint( expr =   400 * model.friend1 +  500 * model.friend2 <= 600)
model.constraints = ConstraintList()
for c in constraints:
    model.constraints.add(
        sum(constraint_coef[c][f] * model.invest_frac[f]
            for f in friends) <= constraint_rhs[c])

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###
print(f"Maximum Profit = ${model.profit():,.2f}")

for f in friends:
    print(f"Batches of {f} = {model.invest_frac[f]():.2f}")

Maximum Profit = $12,000.00
Batches of Friend1 = 0.67
Batches of Friend2 = 0.67

profit_rate

{'Friend1': 9000, 'Friend2': 9000}

Self Assessment: A Holiday Factory

from pyomo.environ import *
import pandas as pd

### PROBLEM DATA ###

# Load Data
num_months = 11
max_duration = 11
months = range(1, num_months + 1)
durations = range(1, max_duration + 1)
rate = [ sum( 20 - i for i in range(month) ) for month in range(1, 12) ]
space_req =  [2000, 2000, 3000, 4000, 6000, 10000, 10000, 10000, 9000, 7000, 5000]

# Parse Dictionaries
rent = dict( zip( durations, rate))
space = dict( zip( months, space_req))

### MODEL CONSTRUCTION ###

# Declaration
model = ConcreteModel(name="HolidayFactory")

# Decision Variables
model.x_sqft = Var(months, durations, domain=NonNegativeReals)

# Objective
model.obj = Objective(expr=sum(rent[d] * model.x_sqft[m, d] for m in months
                               for d in durations))

# Constraints
model.space_ct = ConstraintList()
for month in months:
    model.space_ct.add(
        sum(model.x_sqft[m, d] for m in months for d in durations
            if m <= month and m + d > month) >= space[month])

model.time_rule_ct = ConstraintList()
for m in months:
    for d in durations:
        if m + d > num_months + 1:
            model.time_rule_ct.add(model.x_sqft[m, d] == 0)

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###
print(f"Total Cost = ${model.obj():,.2f}")

print("\nHere are the amounts to lease by month and duration:")
for m in months:
    for d in durations:
        if model.x_sqft[m, d]() > 0:
            print(f"Lease {model.x_sqft[m, d]():.0f} sq ft in month {m:d} for {d:d} months")

Total Cost = $1,125,000.00

Here are the amounts to lease by month and duration:
Lease 2000 sq ft in month 1 for 11 months
Lease 1000 sq ft in month 3 for 9 months
Lease 1000 sq ft in month 4 for 8 months
Lease 1000 sq ft in month 5 for 6 months
Lease 1000 sq ft in month 5 for 7 months
Lease 1000 sq ft in month 6 for 3 months
Lease 2000 sq ft in month 6 for 4 months
Lease 1000 sq ft in month 6 for 5 months

Self Assessment: Supply and Demand Problem

Let $F$ be the set of factories and let $C$ be the set of customers.

Decision Variables: let $x_{f,c}$ be the number of units shipped from factory $f \in F$ to customer $c \in C$

Constants:

• $q_{f,c}$ is the shipping cost per unit between factory $f \in F$ and customer $c \in C$

• $d_c$ is the number of units demanded by customer $c \in C$

• $s_f$ is the number of units supplied by factory $f \in F$

Objective Function: minimize $Cost = \displaystyle \sum_{f \in F} \sum_{c \in C} q_{f,c} x_{f,c}$

Constraints:

• Supply: $\displaystyle \sum_{c \in C} x_{f,c} = s_f, \mbox{ for each } f \in F$

• Demand: $\displaystyle \sum_{f \in F} x_{f,c} = d_c, \mbox{ for each } c \in C$

• Nonnegativity: $x_{f,c} \geq 0$ for each $f \in F, c \in C$

from pyomo.environ import *
import pandas as pd

### PROBLEM DATA

# Load Data
factories = ['factory1', 'factory2']
customers = ['cust1', 'cust2', 'cust3']
usc = [[600, 800, 700], [400, 900, 600]]
sup = [400, 500]
dem = [300, 200, 400]

# Parse into Dictionaries
supply = dict(zip(factories, sup))
demand = dict(zip(customers, dem))
unit_ship_cost = { f:{ c:usc[i][j] for j,c in enumerate(customers)} for i,f in enumerate(factories)}

### MODEL CONSTRUCTION ###

# Declaration
model = ConcreteModel()

# Decision Variables
model.transp = Var(factories, customers, domain=NonNegativeReals)

# Objective
model.total_cost = Objective(expr=sum(unit_ship_cost[f][c] * model.transp[f, c]
                                      for f in factories for c in customers),
                             sense=minimize)

# Constraints
model.supply_ct = ConstraintList()
for f in factories:
    model.supply_ct.add(
        sum(model.transp[f, c] for c in customers) == supply[f])

model.demand_ct = ConstraintList()
for c in customers:
    model.demand_ct.add(
        sum(model.transp[f, c] for f in factories) == demand[c])

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###
print(f"Minimum Total Cost = ${model.total_cost():,.2f}")

# dataframes are displayed nicely in Jupyter
dvars = pd.DataFrame([[model.transp[f, c]() for c in customers]
                      for f in factories],
                     index=factories,
                     columns=customers)
print("Number to ship from each factory to each customer:")
dvars

Minimum Total Cost = $540,000.00
Number to ship from each factory to each customer:

Self Assessment: Positive Shadow Price

Answer: True

Self Assessment: Allowable Range (Objective Coef)

Answer: True

Self Assessment: Changing Parameters

Answer: False

Self Assessment: Graphical Exploration of Sensitivity

(a) Optimal $Z=22$ occurs when $x=6$ and $y=2$.

(b) Shadow price $=23-22=1$. The new optimal is $Z=23$ with coordinates $x=4$ and $y=3$.

(c) The allowable range for resource 2 is $10 \leq b_2 \leq 15$.

The lower bound is 10, as shown here.

The upper bound is 15, as shown here.

(d) The allowable range for the unit profit of activity 2 is $4 \leq c_2 \leq 6$.

The lower bound is 4, as shown here.

The upper bound is 6, as shown here.

Self-Assessment: Solve and Perform Sensitivity

(c) See the code in the following cells and the associated output.

(b) From the column labeled "Marginal" in the top table of the GLPK sensitivity report below, the shadow prices are 0.667 for resource 1 and 1 for resource 2.

from pyomo.environ import *

### PROBLEM DATA ###

# load data
var_names = ['x1', 'x2', 'x3', 'x4']
con_names = ['con1', 'con2']
obj_coef = [5, 4, -1, 3]
con_coef = [ [3, 2, -3, 1], [3, 3, 1, 3] ]
con_rhs = [24, 36]

# parse into dictionaries
obj_coef_dict = dict( zip( var_names, obj_coef) )
con_coef_dict = { c: {v: con_coef[i][j] for j,v in enumerate(var_names) } for i,c in enumerate(con_names)}
con_rhs_dict = dict( zip( con_names, con_rhs))

### MODEL CONSTRUCTION ###

# declaration
model = ConcreteModel(name = "Generic")

# Decision Variables
model.x = Var( var_names, domain = NonNegativeReals)

# Objective
model.obj = Objective( expr = sum( obj_coef_dict[v] * model.x[v] for v in var_names), sense = maximize)

# Constraints
model.cts = ConstraintList()
for c in con_names:
    model.cts.add( sum( con_coef_dict[c][v] * model.x[v] for v in var_names) <= con_rhs_dict[c] )


### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###
print(f"Maximum Z = {model.obj()}")
for v in var_names:
    print(f"{v} = {model.x[v]()}")

Maximum Z = 52.0
x1 = 11.0
x2 = 0.0
x3 = 3.0
x4 = 0.0

# write the model to a sensitivity report
model.write('model.lp', io_options={'symbolic_solver_labels': True})
!glpsol -m model.lp --lp --ranges sensit.sen

GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 -m model.lp --lp --ranges sensit.sen
Reading problem data from 'model.lp'...
2 rows, 4 columns, 8 non-zeros
31 lines were read
GLPK Simplex Optimizer 5.0
2 rows, 4 columns, 8 non-zeros
Preprocessing...
2 rows, 4 columns, 8 non-zeros
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  3.000e+00  ratio =  3.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 2
*     0: obj =  -0.000000000e+00 inf =   0.000e+00 (3)
*     2: obj =   5.200000000e+01 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.0 Mb (39693 bytes)
Write sensitivity analysis report to 'sensit.sen'...

# widen browser and/or close TOC to see sensitivity report
import numpy as np
np.set_printoptions(linewidth=110)
f = open('sensit.sen', 'r')
file_contents = f.read()
print(file_contents)
f.close()

GLPK 5.0  - SENSITIVITY ANALYSIS REPORT                                                                         Page   1

Problem:    
Objective:  obj = 52 (MAXimum)

   No. Row name     St      Activity         Slack   Lower bound       Activity      Obj coef  Obj value at Limiting
                                          Marginal   Upper bound          range         range   break point variable
------ ------------ -- ------------- ------------- -------------  ------------- ------------- ------------- ------------
     1 c_u_cts(1)_  NU      24.00000        .               -Inf     -108.00000       -.66667     -36.00000 x(x1)
                                            .66667      24.00000       36.00000          +Inf      60.00000 x(x3)

     2 c_u_cts(2)_  NU      36.00000        .               -Inf       24.00000      -1.00000      40.00000 x(x3)
                                           1.00000      36.00000           +Inf          +Inf          +Inf

GLPK 5.0  - SENSITIVITY ANALYSIS REPORT                                                                         Page   2

Problem:    
Objective:  obj = 52 (MAXimum)

   No. Column name  St      Activity      Obj coef   Lower bound       Activity      Obj coef  Obj value at Limiting
                                          Marginal   Upper bound          range         range   break point variable
------ ------------ -- ------------- ------------- -------------  ------------- ------------- ------------- ------------
     1 x(x1)        BS      11.00000       5.00000        .              .            4.63636      48.00000 x(x2)
                                            .               +Inf       11.00000          +Inf          +Inf

     2 x(x2)        NL        .            4.00000        .                -Inf          -Inf          +Inf
                                           -.33333          +Inf       12.00000       4.33333      48.00000 x(x1)

     3 x(x3)        BS       3.00000      -1.00000        .            -3.60000      -2.33333      48.00000 x(x4)
                                            .               +Inf       36.00000       1.66667      60.00000 c_u_cts(1)_

     4 x(x4)        NL        .            3.00000        .                -Inf          -Inf          +Inf
                                           -.66667          +Inf        6.00000       3.66667      48.00000 x(x3)

End of report

Shadow Prices

• Resource 1: 0.66667

• Resource 2: 1.00000

Allowable Range for Right Hand Side of Constraints

• Resource 1: -108 to 36

• Resource 2: 24 to $\infty$

Allowable Range Objective Function Coefficients

• Coefficient 1: 4.63636 to $\infty$

• Coefficient 2: $-\infty$ to 4.3333

• Coefficient 3: -2.33333 to 1.66667

• Coefficient 4: $-\infty$ to 3.66667

Self-Assessment: Formulate, Solve, and Perform Sensitivity #1

(a) See the code and output for the cell below. The maximum profit is $3500, obtained when 2000 toys and 1000 subassemlies are produced per day.

(f) From the "Obj coef range" in the bottom table of the GLPK sensitivty report, the allowable range of the unit profit for toys is $2.50 to $5 whereas that for subassemblies is -$3 to -$1.50.

from pyomo.environ import *

### PROBLEM DATA ###

# load data
parts = ['toys','subs']
unit_profit = [3,-2.5]
con_names = ['subA','subB']
con_rhs = [3000, 1000]
con_coef = [ [2,-1], [1,-1] ]

# parse dictionaries
unit_profit_dict = dict( zip( parts, unit_profit))
con_rhs_dict = dict( zip( con_names, con_rhs) )
con_coef_dict = { c: {p:con_coef[i][j] for j,p in enumerate(parts)} for i,c in enumerate(con_names) }

### MODEL CONSTRUCTION ###

# Declaration
model = ConcreteModel(name = "TannerCo")

# Decision Variables
model.x = Var( parts, domain = NonNegativeReals)

# Objective
model.obj = Objective( expr = sum( unit_profit_dict[p] * model.x[p] for p in parts), sense = maximize )

# Constraints
model.cts = ConstraintList()
for c in con_names:
    model.cts.add( sum( con_coef_dict[c][p] * model.x[p] for p in parts) <= con_rhs_dict[c] )

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###
print(f"Maximum Profit = ${model.obj():,.2f}")
for p in parts:
    print(f"Daily product rate of {p} = {model.x[p]():.0f}")

Maximum Profit = $3,500.00
Daily product rate of toys = 2000
Daily product rate of subs = 1000

# write the model to a sensitivity report
model.write('model.lp', io_options={'symbolic_solver_labels': True})
!glpsol -m model.lp --lp --ranges sensit.sen

GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 -m model.lp --lp --ranges sensit.sen
Reading problem data from 'model.lp'...
2 rows, 2 columns, 4 non-zeros
23 lines were read
GLPK Simplex Optimizer 5.0
2 rows, 2 columns, 4 non-zeros
Preprocessing...
2 rows, 2 columns, 4 non-zeros
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  2.000e+00  ratio =  2.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 2
*     0: obj =  -0.000000000e+00 inf =   0.000e+00 (1)
*     2: obj =   3.500000000e+03 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.0 Mb (39693 bytes)
Write sensitivity analysis report to 'sensit.sen'...

# widen browser and/or close TOC to see sensitivity report
import numpy as np
np.set_printoptions(linewidth=110)
f = open('sensit.sen', 'r')
file_contents = f.read()
print(file_contents)
f.close()

GLPK 4.65 - SENSITIVITY ANALYSIS REPORT                                                                         Page   1

Problem:    
Objective:  obj = 3500 (MAXimum)

   No. Row name     St      Activity         Slack   Lower bound       Activity      Obj coef  Obj value at Limiting
                                          Marginal   Upper bound          range         range   break point variable
------ ------------ -- ------------- ------------- -------------  ------------- ------------- ------------- ------------
     1 c_u_cts(1)_  NU    3000.00000        .               -Inf     2000.00000       -.50000    3000.00000 x(subs)
                                            .50000    3000.00000           +Inf          +Inf          +Inf

     2 c_u_cts(2)_  NU    1000.00000        .               -Inf           -Inf      -2.00000          -Inf
                                           2.00000    1000.00000     1500.00000          +Inf    4500.00000 x(subs)

     3 c_e_ONE_VAR_CONSTANT
                    NS       1.00000        .            1.00000         .               -Inf    3500.00000 ONE_VAR_CONSTANT
                                            .            1.00000           +Inf          +Inf    3500.00000

GLPK 4.65 - SENSITIVITY ANALYSIS REPORT                                                                         Page   2

Problem:    
Objective:  obj = 3500 (MAXimum)

   No. Column name  St      Activity      Obj coef   Lower bound       Activity      Obj coef  Obj value at Limiting
                                          Marginal   Upper bound          range         range   break point variable
------ ------------ -- ------------- ------------- -------------  ------------- ------------- ------------- ------------
     1 x(subs)      BS    1000.00000      -2.50000        .         -1000.00000      -3.00000    3000.00000 c_u_cts(1)_
                                            .               +Inf           +Inf      -1.50000    4500.00000 c_u_cts(2)_

     2 x(toys)      BS    2000.00000       3.00000        .          1000.00000       2.50000    2500.00000 c_u_cts(1)_
                                            .               +Inf           +Inf       5.00000    7500.00000 c_u_cts(2)_

     3 ONE_VAR_CONSTANT
                    BS       1.00000        .             .             1.00000          -Inf          -Inf
                                            .               +Inf        1.00000          +Inf          +Inf

End of report

Self-Assessment: Formulate, Solve, and Perform Sensitivity #2

(a) This is pretty simple. Add an extra constraint that $x_1 \leq 2500$. The maximum profit is still $3500, obtained when 2000 toys and 1000 subassemlies are produced per day. See the code cell below.

See cells befow for parts b and f.

from pyomo.environ import *

### PROBLEM DATA ###

# load data
parts = ['toys','subs']
unit_profit = [3,-2.5]
con_names = ['subA','subB','max_toys']
con_rhs = [3000, 1000, 2500]
con_coef = [ [2,-1], [1,-1], [1,0] ]

# parse dictionaries
unit_profit_dict = dict( zip( parts, unit_profit))
con_rhs_dict = dict( zip( con_names, con_rhs) )
con_coef_dict = { c: {p:con_coef[i][j] for j,p in enumerate(parts)} for i,c in enumerate(con_names) }

### MODEL CONSTRUCTION ###

# Declaration
model = ConcreteModel(name = "TannerCo")

# Decision Variables
model.x = Var( parts, domain = NonNegativeReals)

# Objective
model.obj = Objective( expr = sum( unit_profit_dict[p] * model.x[p] for p in parts), sense = maximize )

# Constraints
model.cts = ConstraintList()
for c in con_names:
    model.cts.add( sum( con_coef_dict[c][p] * model.x[p] for p in parts) <= con_rhs_dict[c] )

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###
print(f"Maximum Profit = ${model.obj():,.2f}")
for p in parts:
    print(f"Daily product rate of {p} = {model.x[p]():.0f}")

Maximum Profit = $3,500.00
Daily product rate of toys = 2000
Daily product rate of subs = 1000

(b) Run the model with RHS of coefficient 1 at 3001 rather than 3000. The shadow price for subassembly A is $0.50, which is the maximum premium that the company should be willing to pay. See the code in the next cell.

from pyomo.environ import *

### PROBLEM DATA ###

# load data
parts = ['toys','subs']
unit_profit = [3,-2.5]
con_names = ['subA','subB','max_toys']
con_rhs = [3001, 1000, 2500]
con_coef = [ [2,-1], [1,-1], [1,0] ]

# parse dictionaries
unit_profit_dict = dict( zip( parts, unit_profit))
con_rhs_dict = dict( zip( con_names, con_rhs) )
con_coef_dict = { c: {p:con_coef[i][j] for j,p in enumerate(parts)} for i,c in enumerate(con_names) }

### MODEL CONSTRUCTION ###

# Declaration
model = ConcreteModel(name = "TannerCo")

# Decision Variables
model.x = Var( parts, domain = NonNegativeReals)

# Objective
model.obj = Objective( expr = sum( unit_profit_dict[p] * model.x[p] for p in parts), sense = maximize )

# Constraints
model.cts = ConstraintList()
for c in con_names:
    model.cts.add( sum( con_coef_dict[c][p] * model.x[p] for p in parts) <= con_rhs_dict[c] )

### SOLUTION ###
solver = SolverFactory('glpk')
solver.solve(model)

### OUTPUT ###
print(f"Maximum Profit = ${model.obj():,.2f}")
for p in parts:
    print(f"Daily product rate of {p} = {model.x[p]():.0f}")

Maximum Profit = $3,500.50
Daily product rate of toys = 2001
Daily product rate of subs = 1001

(f) As shown in the sensitivity report, the shadow price is $0.50 for subassembly A and $2 for subassembly B. According to the activity range, the allowable range for the right-hand side of the subassembly A constraint is 2,000 to 3,500. The allowable range for the right-hand side of the subassembly B constraint is 500 to 1,500.

# write the model to a sensitivity report
model.write('model.lp', io_options={'symbolic_solver_labels': True})
!glpsol -m model.lp --lp --ranges sensit.sen

# widen browser and/or close TOC to see sensitivity report
import numpy as np
np.set_printoptions(linewidth=110)
f = open('sensit.sen', 'r')
file_contents = f.read()
print(file_contents)
f.close()

GLPSOL--GLPK LP/MIP Solver 5.0
Parameter(s) specified in the command line:
 -m model.lp --lp --ranges sensit.sen
Reading problem data from 'model.lp'...
3 rows, 2 columns, 5 non-zeros
27 lines were read
GLPK Simplex Optimizer 5.0
3 rows, 2 columns, 5 non-zeros
Preprocessing...
2 rows, 2 columns, 4 non-zeros
Scaling...
 A: min|aij| =  1.000e+00  max|aij| =  2.000e+00  ratio =  2.000e+00
Problem data seem to be well scaled
Constructing initial basis...
Size of triangular part is 2
*     0: obj =  -0.000000000e+00 inf =   0.000e+00 (1)
*     2: obj =   3.500500000e+03 inf =   0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used:   0.0 secs
Memory used: 0.0 Mb (40412 bytes)
Write sensitivity analysis report to 'sensit.sen'...
GLPK 5.0  - SENSITIVITY ANALYSIS REPORT                                                                         Page   1

Problem:    
Objective:  obj = 3500.5 (MAXimum)

   No. Row name     St      Activity         Slack   Lower bound       Activity      Obj coef  Obj value at Limiting
                                          Marginal   Upper bound          range         range   break point variable
------ ------------ -- ------------- ------------- -------------  ------------- ------------- ------------- ------------
     1 c_u_cts(1)_  NU    3001.00000        .               -Inf     2000.00000       -.50000    3000.00000 x(subs)
                                            .50000    3001.00000     3500.00000          +Inf    3750.00000 c_u_cts(3)_

     2 c_u_cts(2)_  NU    1000.00000        .               -Inf      501.00000      -2.00000    2502.50000 c_u_cts(3)_
                                           2.00000    1000.00000     1500.50000          +Inf    4501.50000 x(subs)

     3 c_u_cts(3)_  BS    2001.00000     499.00000          -Inf     1000.00000       -.50000    2500.00000 c_u_cts(1)_
                                            .         2500.00000           +Inf       2.00000    7502.50000 c_u_cts(2)_

GLPK 5.0  - SENSITIVITY ANALYSIS REPORT                                                                         Page   2

Problem:    
Objective:  obj = 3500.5 (MAXimum)

   No. Column name  St      Activity      Obj coef   Lower bound       Activity      Obj coef  Obj value at Limiting
                                          Marginal   Upper bound          range         range   break point variable
------ ------------ -- ------------- ------------- -------------  ------------- ------------- ------------- ------------
     1 x(toys)      BS    2001.00000       3.00000        .          1000.00000       2.50000    2500.00000 c_u_cts(1)_
                                            .               +Inf     2500.00000       5.00000    7502.50000 c_u_cts(2)_

     2 x(subs)      BS    1001.00000      -2.50000        .         -1000.00000      -3.00000    3000.00000 c_u_cts(1)_
                                            .               +Inf     1999.00000      -1.50000    4501.50000 c_u_cts(2)_

End of report