The LOCAL JOB SHOP problem from Hillier and Lieberman, Section 11.3, can also be solved "brute-force" by using a solver for generic (convex) nonlinear programs such as Ipopt.

from pyomo.environ import *
from pyomo.opt import *
opt = solvers.SolverFactory("ipopt")

model = ConcreteModel()

To simplify the code, we use 5 decision variables, where the first and the last are the spring level of employment which is constrained to a single number anyway. The advantage is that we can formulate the objective function without explicit special cases. $S$ is our index set and $r$ is the vector containing the minimum required employment levels.

S= range(5)
r = [255, 220, 240, 200, 255]

Instead of introducing explicit constraints, we can put the constraints on the employment levels as bounds on the range of the variable directly into the definition of the decision variable.

def f(model, i):
    return (r[i],255)

model.x = Var(S, bounds=f)

Then we only need to define the objective function as the sum of all partial costs.

model.c = Objective(expr = sum(200*(model.x[i]-model.x[i+1])**2 + 
                               2000*(model.x[i]-r[i]) for i in range(4)),
                    sense=minimize)

results = opt.solve(model)

The resulting optimal value are, of course, the same as result from the dynamic programming analysis.

model.x.get_values()

{0: 255.0,
 1: 247.49999999988543,
 2: 244.99999999998957,
 3: 247.49999999986818,
 4: 255.0}

model.c.expr()

185000.0