<jupyter_start><jupyter_text>![image.png](attachment:image.png)<jupyter_code>from pyomo.environ import *
from math import pi<jupyter_output><empty_output><jupyter_text>Concrete Pyomo Model with Explicit Variables<jupyter_code>m = ConcreteModel()

m.x_1 = Var(within=NonNegativeReals)
m.x_2 = Var(within=NonNegativeReals)

m.obj = Objective(expr= m.x_1 + 2*m.x_2, sense=minimize )

m.c1 = Constraint(expr= 3*m.x_1 + 4*m.x_2 >= 1)
m.c2 = Constraint(expr= 2*m.x_1 + 5*m.x_2 >= 2)

solver = SolverFactory('ipopt')
status = solver.solve(m)

print("Status = %s" % status.solver.termination_condition)

print("%s = %f" % (m.x_1, value(m.x_1)))
print("%s = %f" % (m.x_2, value(m.x_2)))
print("Objective = %f" % value(m.obj))<jupyter_output>Status = optimal
x_1 = 0.000000
x_2 = 0.400000
Objective = 0.800000<jupyter_text>Concrete Pyomo Model with Indexed Variables<jupyter_code>model = ConcreteModel()
model.x = Var([1,2], within=NonNegativeReals)
model.obj = Objective(expr=model.x[1] + 2*model.x[2])
model.con1 = Constraint(expr=3*model.x[1] + 4*model.x[2]>=1)
model.con2 = Constraint(expr=2*model.x[1] + 5*model.x[2]>=2)

solver = SolverFactory('ipopt')
status = solver.solve(model)

print("Status = %s" % status.solver.termination_condition)

print("%s = %f" % (model.x[1], value(model.x[1])))
print("%s = %f" % (model.x[2], value(model.x[2])))
print("Objective = %f" % value(model.obj))
# Modificando solo el nombre de las variables 
model = ConcreteModel()
model.x = Var([1,3], within=NonNegativeReals)

model.obj = Objective(expr=model.x[1] + 2*model.x[3])

model.con1 = Constraint(expr=3*model.x[1] + 4*model.x[3]>=1)
model.con2 = Constraint(expr=2*model.x[1] + 5*model.x[3]>=2)

solver = SolverFactory('ipopt')
status = solver.solve(model)

print("Status = %s" % status.solver.termination_condition)

print("%s = %f" % (model.x[1], value(model.x[1])))
print("%s = %f" % (model.x[3], value(model.x[3])))
print("Objective = %f" % value(model.obj))<jupyter_output><empty_output><jupyter_text>Concrete Pyomo Model with External Data ![image.png](attachment:image.png)<jupyter_code>N = [1,2] # Nombre de las variables

c = {1:1, 2:2} # Costos para la funcion de costo 
# Ej. si fuera 5x_1 + 6x_2  tendriamos c = {1:5, 2:6}

a = {(1,1):3, (2,1):4, (1,2):2, (2,2):5} # Parametros de suget to s.t 
# Es una matriz

b = {1:1, 2:2} #  bounds  o cotas

model = ConcreteModel()
model.x = Var(N, within=NonNegativeReals)

model.obj = Objective(expr=sum(c[i]*model.x[i] for i in N))

model.con1 = Constraint(expr=sum(a[i,1]*model.x[i] for i in N) >= b[1])
model.con2 = Constraint(expr=sum(a[i,2]*model.x[i] for i in N) >= b[2])

solver = SolverFactory('ipopt')
status = solver.solve(model)

print("Status = %s" % status.solver.termination_condition)

for i in N:
    print("%s = %f" % (model.x[i], value(model.x[i])))

print("Objective = %f" % value(model.obj))<jupyter_output>Status = optimal
x[1] = 0.000000
x[2] = 0.400000
Objective = 0.800000<jupyter_text>Concrete Pyomo Model with Constraint Rules<jupyter_code>N = [1,2]
M = [1,2]
c = {1:1, 2:2}
a = {(1,1):3, (2,1):4, (1,2):2, (2,2):5}
b = {1:1, 2:2}

model = ConcreteModel()

model.x = Var(N, within=NonNegativeReals)
model.obj = Objective(expr=sum(c[i]*model.x[i] for i in N))

def con_rule(model, m):
    return sum(a[i,m]*model.x[i] for i in N) >= b[m]

model.con = Constraint(M, rule=con_rule)


solver = SolverFactory('ipopt')
status = solver.solve(model)

print("Status = %s" % status.solver.termination_condition)

for i in N:
    print("%s = %f" % (model.x[i], value(model.x[i])))

print("Objective = %f" % value(model.obj))<jupyter_output>Status = optimal
x[1] = 0.000000
x[2] = 0.400000
Objective = 0.800000<jupyter_text>Concrete Pyomo Model with Abstract Component Declarations<jupyter_code>model = ConcreteModel()

def N_rule(model):
    return [1,2]

model.N = Set(rule=N_rule)
model.M = Set(initialize=[1,2])

model.c = Param(model.N, initialize={1:1, 2:2})
model.a = Param(model.N, model.M,initialize={(1,1):3, (2,1):4, (1,2):2, (2,2):5})
model.b = Param(model.M, initialize={1:1, 2:2})

model.x = Var(model.N, within=NonNegativeReals)

# Objective
def obj_rule(model):
    return sum(model.c[i]*model.x[i] for i in model.N)
model.obj = Objective(rule=obj_rule)

# Constraint
def con_rule(model, m):
    return sum(model.a[i,m]*model.x[i] for i in model.N)  >= model.b[m]

model.con = Constraint(model.M, rule=con_rule)



solver = SolverFactory('ipopt')
status = solver.solve(model)

print("Status = %s" % status.solver.termination_condition)

for i in N:
    print("%s = %f" % (model.x[i], value(model.x[i])))

print("Objective = %f" % value(model.obj))<jupyter_output>Status = optimal
x[1] = 0.000000
x[2] = 0.400000
Objective = 0.800000<jupyter_text>Using a Function to Construct a Concrete Pyomo Model<jupyter_code>def create_model(N=[], M=[], c={}, a={}, b={}):
    model = ConcreteModel()
    model.x = Var(N, within=NonNegativeReals)
    model.obj = Objective(expr=sum(c[i]*model.x[i] for i in N))

    def con_rule(model, m):
        return sum(a[i,m]*model.x[i] for i in N) >= b[m]
    model.con = Constraint(M, rule=con_rule)
    return model

model = create_model(N = [1,2], M = [1,2], c = {1:1, 2:2},
        a = {(1,1):3, (2,1):4, (1,2):2, (2,2):5},
        b = {1:1, 2:2})

solver = SolverFactory('ipopt')
status = solver.solve(model)

print("Status = %s" % status.solver.termination_condition)

for i in N:
    print("%s = %f" % (model.x[i], value(model.x[i])))

print("Objective = %f" % value(model.obj))<jupyter_output>Status = optimal
x[1] = 0.000000
x[2] = 0.400000
Objective = 0.800000<jupyter_text>Abstract Pyomo Model<jupyter_code>model = AbstractModel()

model.N = Set()
model.M = Set()

model.c = Param(model.N)
model.a = Param(model.N, model.M)
model.b = Param(model.M)

model.x = Var(model.N, within=NonNegativeReals)

def obj_rule(model):
    return sum(model.c[i]*model.x[i] for i in model.N)
model.obj = Objective(rule=obj_rule)

def con_rule(model, m):
    return sum(model.a[i,m]*model.x[i] for i in model.N) >= model.b[m]
model.con = Constraint(model.M, rule=con_rule)
%%writefile info.dat

set N := 1 2 ;
    
set M := 1 2 ;
    
param c :=
1 1
2 2;

param a :=
1 1 3
2 1 4
1 2 2
2 2 5;
param b :=
1 1
2 2 ;
# %load info.dat

set N := 1 2 ;
    
set M := 1 2 ;
    
param c :=
1 1
2 2;

param a :=
1 1 3
2 1 4
1 2 2
2 2 5;
param b :=
1 1
2 2 ;
instance = model.create_instance('info.dat')

#instance.pprint()
opt = SolverFactory("ipopt")
results = opt.solve(instance)
results.write()

solver = SolverFactory('glpk')

instance = model.create_instance('info.dat')
#results = solver.solve(instance,tee=True)
results = solver.solve(instance)
results.write()
instance.solutions.load_from(results)

for v in instance.component_objects(Var, active=True):
    print ("Variable",v)
    varobject = getattr(instance, str(v))
    for index in varobject:
        print ("   ",index, varobject[index].value)
instance.solutions.load_from(results)

for v in instance.component_objects(Var, active=True):
    print ("Variable",v)
    varobject = getattr(instance, str(v))
    for index in varobject:
        print ("   ",index, varobject[index].value)<jupyter_output>Variable x
    1 0.0
    2 0.4<jupyter_text>![image.png](attachment:image.png) http://cvxopt.org/examples/tutorial/lp.html ![image.png](attachment:image.png)<jupyter_code>from cvxopt import matrix, solvers

A = matrix([ [-1.0, -1.0, 0.0, 1.0], [1.0, -1.0, -1.0, -2.0] ])

b = matrix([ 1.0, -2.0, 0.0, 4.0 ])

c = matrix([ 2.0, 1.0 ])

sol=solvers.lp(c,A,b)
print(sol['x'])<jupyter_output>[ 5.00e-01]
[ 1.50e+00]<jupyter_text>![image.png](attachment:image.png)<jupyter_code>A = matrix([ [-3.0, -2.0, -1.0, 0.0],
             [-4.0, -5.0, 0.0, -1.0] ])

b = matrix([ -1.0, -2.0, 0.0, 0 ])

c = matrix([ 1.0, 2.0 ])

sol=solvers.lp(c,A,b)

print(sol['x'])
# Abstract Knapsack Problem
#

from pyomo.environ import *

model = AbstractModel()

model.ITEMS = Set()

model.v = Param(model.ITEMS, within=PositiveReals)

model.w = Param(model.ITEMS, within=PositiveReals)

model.limit = Param(within=PositiveReals)

model.x = Var(model.ITEMS, within=Binary)

def value_rule(model):
    return sum(model.v[i]*model.x[i] for i in model.ITEMS)
model.value = Objective(sense=maximize, rule=value_rule)

def weight_rule(model):
    return sum(model.w[i]*model.x[i] for i in model.ITEMS) <= model.limit
model.weight = Constraint(rule=weight_rule)



data = {
    'ITEMS': {None:('hammer','wrench','screwdriver','towel')},
    'v': {'hammer': 8,
           'wrench': 3,
           'screwdriver': 6,
           'towel': 11,
       },
    'w': { 'hammer': 5,
           'wrench': 7,
           'screwdriver': 4,
           'towel': 3,
       },
    'limit': {None:14},
}
inst = model.create_instance(data={None:data})
inst.pprint()

solver = SolverFactory('glpk')
results = solver.solve(inst)
results.write()
instance.solutions.load_from(results)

for v in instance.component_objects(Var, active=True):
    print ("Variable",v)
    varobject = getattr(instance, str(v))
    for index in varobject:
        print ("   ",index, varobject[index].value)
%matplotlib inline
# Example 1 from http://www.mathworks.com/help/matlab/ref/pdepe.html

from pyomo.environ import *
from pyomo.dae import *
from pyomo.opt import SolverFactory
from pyomo.dae.plugins.finitedifference import Finite_Difference_Transformation
from pyomo.dae.plugins.colloc import Collocation_Discretization_Transformation
import math

m = ConcreteModel()
m.t = ContinuousSet(bounds=(0,2))
m.x = ContinuousSet(bounds=(0,1))
m.u = StateVar(m.x,m.t)

m.dudx = DerivativeVar(m.u,wrt=m.x)
m.dudx2 = DerivativeVar(m.u,wrt=(m.x,m.x))
m.dudt = DerivativeVar(m.u,wrt=m.t)

def _pde(m,i,j):
    if i == 0 or i == 1 or j == 0 :
        return Constraint.Skip
    return math.pi**2*m.dudt[i,j] == m.dudx2[i,j]
m.pde = Constraint(m.x,m.t,rule=_pde)

def _initcon(m,i):
    if i == 0 or i == 1:
        return Constraint.Skip
    return m.u[i,0] == sin(math.pi*i)
m.initcon = Constraint(m.x,rule=_initcon)

def _lowerbound(m,j):
    return m.u[0,j] == 0
m.lowerbound = Constraint(m.t,rule=_lowerbound)

def _upperbound(m,j):
    return math.pi*exp(-j)+m.dudx[1,j] == 0
m.upperbound = Constraint(m.t,rule=_upperbound)

m.obj = Objective(expr=1)

# Discretize using Finite Difference Method
discretize = Finite_Difference_Transformation()
disc = discretize.apply(m,nfe=25,wrt=m.x,scheme='BACKWARD')
disc = discretize.apply(disc,nfe=20,wrt=m.t,scheme='BACKWARD',clonemodel=False)

# Discretize using Orthogonal Collocation
#discretize2 = Collocation_Discretization_Transformation()
#disc = discretize2.apply(disc,nfe=10,ncp=3,wrt=m.x,clonemodel=False)
#disc = discretize2.apply(disc,nfe=20,ncp=3,wrt=m.t,clonemodel=False)


solver='ipopt'
opt=SolverFactory(solver)

results = opt.solve(disc,tee=True)
disc.load(results)

#disc.u.pprint()

x = []
t = []
u = []

for i in sorted(disc.x):
    temp=[]
    tempx = []
    for j in sorted(disc.t):
        tempx.append(i)
        temp.append(value(disc.u[i,j]))
    x.append(tempx)
    t.append(sorted(disc.t))
    u.append(temp)


import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(1,1,1,projection='3d')
ax.set_xlabel('Distance x')
ax.set_ylabel('Time t')
ax.set_title('Numerical Solution Using Backward Difference Method')
p = ax.plot_wireframe(x,t,u,rstride=1,cstride=1)
fig.show()
from pyomo.environ import *

# Creates a list of the Ingredients
Ingredients = ['CHICKEN', 'BEEF', 'MUTTON', 'RICE', 'WHEAT', 'GEL']

# A dictionary of the costs of each of the Ingredients is created
costs = {'CHICKEN': 0.013,
         'BEEF': 0.008,
         'MUTTON': 0.010,
         'RICE': 0.002,
         'WHEAT': 0.005,
         'GEL': 0.001}

# A dictionary of the protein percent in each of the Ingredients is created
proteinPercent = {'CHICKEN': 0.100,
                  'BEEF': 0.200,
                  'MUTTON': 0.150,
                  'RICE': 0.000,
                  'WHEAT': 0.040,
                  'GEL': 0.000}

# A dictionary of the fat percent in each of the Ingredients is created
fatPercent = {'CHICKEN': 0.080,
              'BEEF': 0.100,
              'MUTTON': 0.110,
              'RICE': 0.010,
              'WHEAT': 0.010,
              'GEL': 0.000}

# A dictionary of the fibre percent in each of the Ingredients is created
fibrePercent = {'CHICKEN': 0.001,
                'BEEF': 0.005,
                'MUTTON': 0.003,
                'RICE': 0.100,
                'WHEAT': 0.150,
                'GEL': 0.000}

# A dictionary of the salt percent in each of the Ingredients is created
saltPercent = {'CHICKEN': 0.002,
               'BEEF': 0.005,
               'MUTTON': 0.007,
               'RICE': 0.002,
               'WHEAT': 0.008,
               'GEL': 0.000}

model = ConcreteModel(name="The Whiskas Problem")

model.ingredient_vars = Var(Ingredients, bounds=(0,None), doc="The amount of each ingredient that is used")

model.obj = Objective(expr=sum(costs[i]*model.ingredient_vars[i] for i in Ingredients), doc="Total Cost of Ingredients per can")

model.c0 = Constraint(expr=sum(model.ingredient_vars[i] for i in Ingredients) == 100, doc="PercentagesSum")
model.c1 = Constraint(expr=sum(proteinPercent[i] * model.ingredient_vars[i] for i in Ingredients) >= 8.0, doc="ProteinRequirement")
model.c2 = Constraint(expr=sum(fatPercent[i] * model.ingredient_vars[i] for i in Ingredients) >= 6.0, doc="FatRequirement")
model.c3 = Constraint(expr=sum(fibrePercent[i] * model.ingredient_vars[i] for i in Ingredients) <= 2.0, doc="FibreRequirement")
model.c4 = Constraint(expr=sum(saltPercent[i] * model.ingredient_vars[i] for i in Ingredients) <= 0.4, doc="SaltRequirement")
#inst = model.create_instance(data={None:data})
#inst.pprint()

solver = SolverFactory('glpk')
results = solver.solve(model)
results.write()
instance.solutions.load_from(results)

for v in instance.component_objects(Var, active=True):
    print ("Variable",v)
    varobject = getattr(instance, str(v))
    for index in varobject:
        print ("   ",index, varobject[index].value)
%matplotlib notebook
# sodacan.py
from pyomo.environ import *
from math import pi

M = ConcreteModel()
M.r = Var(bounds=(0,None))
M.h = Var(bounds=(0,None))
M.o = Objective(expr=\
        2*pi*M.r*(M.r + M.h))
M.c = Constraint(expr=\
        pi*M.h*M.r**2 == 355)


solver = SolverFactory('ipopt')
status = solver.solve(M)
status.write()

print("Status = %s" % status.solver.termination_condition)
print("%s = %f" % (M.r, value(M.r)))
print("%s = %f" % (M.h, value(M.h)))

print("Objective = %f" % value(M.o))
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as plt
import numpy as np
from math import pi, sqrt

fig = plt.figure()
ax = fig.gca(projection='3d')
R_ = np.arange(0.25, 10, 0.25)
H_ = np.arange(0.25, 10, 0.25)
R, H = np.meshgrid(R_, H_)
Z = 2*pi*R*(R+H)
surf = ax.plot_surface(R, H, Z, rstride=1, cstride=1, cmap=cm.hot,
                       linewidth=0, antialiased=False)
ax.set_xlabel("r")
ax.set_ylabel("h")
ax.set_zlim(0, 1200)

ax.zaxis.set_major_locator(LinearLocator(10))
#ax.zaxis.set_major_formatter(FormatStrFormatter(' %.02f'))

#fig.colorbar(surf, shrink=0.5, aspect=5)

H_ = 355/(pi*R_*R_)
valid = np.where(H_ < 10.1)
Z_ = R_+H_
Z_ = 2*pi*R_*Z_
ax.plot(R_[valid], H_[valid], Z_[valid], label='parametric curve')
ax.scatter3D(3.837215,7.674430,0,s=20,c="green")
ax.scatter3D?
plt.show()
model = ConcreteModel()

model.x = Var( initialize=-1.2, bounds=(-2, 2) )
model.y = Var( initialize= 1.0, bounds=(-2, 2) )

model.obj = Objective(
expr= (1-model.x)**2 + 100*(model.y-model.x**2)**2,
sense= minimize )
solver = SolverFactory('ipopt')
status = solver.solve(model)
status.write()

print("Status = %s" % status.solver.termination_condition)
print("%s = %f" % (model.x, value(model.x)))
print("%s = %f" % (model.y, value(model.y)))

print("Objective = %f" % value(model.obj))
print("\nDisplaying Solution\n" + '-'*60)
pyomo_postprocess(None, model, status)
## Display of the output ##
# Display x.l, x.m ;
def pyomo_postprocess(options=None, instance=None, results=None):
    model.x.display()
model.x.display()
model.y.display()<jupyter_output>y : Size=1, Index=None, Domain=Reals
    Key  : Lower : Value              : Upper : Fixed : Stale
    None :    -2 : 0.9999999949529991 :     2 : False : False