using JuMP, Ipopt

m = Model(solver = IpoptSolver(print_level=0))
@variable(m, 25 <= r <= 60)
@variable(m, h >= 0)
@variable(m, c >= 50)

@constraint(m, r + h + c == 100)

@NLconstraint(m, tensile,  12.5 - 0.1h - 0.001h^2 >= 12)
@NLconstraint(m, elasticity, 17 + 0.35r - 0.04h - 0.002r^2 >= 16)
@NLconstraint(m, hard_lo, 25 <= 34 + 0.1r + 0.06h - 0.3c + 0.01r*h + 0.005h^2 + 0.001*c^1.95)
@NLconstraint(m, hard_hi, 34 + 0.1r + 0.06h - 0.3c + 0.01r*h + 0.005h^2 + 0.001*c^1.95 <= 35)

@objective(m, Min, 0.04r + 0.01h + 0.07c )
solve(m)

# rubber, oil, and carbon
getvalue([r h c])

getdual(tensile)

getdual(elasticity)

getdual(hard_lo)

getdual(hard_hi)

using JuMP, Ipopt

ph = 0.05
pr = 0.1
pc = 0.02

m = Model(solver = IpoptSolver(print_level=0))
@variable(m, 4.77 >= h >= 0)
@variable(m, c >= 50)

@constraint(m, 50 <= h+c <= 75)

@NLconstraint(m, elasticity, 32 + 0.05c - 0.002c^2 + 0.01h - 0.004c*h - 0.002h^2 >= 16)
@NLconstraint(m, hard_lo, 25 <= 44 + 0.96h - 0.4c - 0.01c*h + 0.005h^2 + 0.001*c^1.95)
@NLconstraint(m, hard_hi, 44 + 0.96h - 0.4c - 0.01c*h + 0.005h^2 + 0.001*c^1.95 <= 35)

@objective(m, Min, (ph-pr)*h + (pc-pr)*c )
solve(m)

# rubber, oil, and carbon
getvalue([100-h-c h c])

getdual(elasticity)

getdual(hard_lo)

getdual(hard_hi)