Advanced topics

Removing the constraint on Beta

We can remove the constrinat on Beta by adding model.__model__.beta.setlb(None) before the model optimization (i.e., model.optimize()). Now the estimated Beta can take any value in Reals. For example, the constraint Beta >= 0 is removed in CQR:

# import packages
from pystoned import CQER
from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
from pystoned import dataset as dataset

# import the GHG example data
data = dataset.load_GHG_abatement_cost(x_select=['HRSN', 'CPNK', 'GHG'], y_select=['VALK'])

# calculate the quantile model
model = CQER.CQR(y=data.y, x=data.x, tau=0.5, z=None, cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
# remove the constraint beta >= 0
model.__model__.beta.setlb(None)
model.optimize(OPT_LOCAL)

# display estimated beta
model.display_beta()

Adding an additional constraint

The extral constraint, e.g., 0 <= Beta <= 1, can be added using

from pyomo.environ import Constraint

def constraint_rule(model, i, j):
    upperbound = [1.0, 1.0, 1.0]
    return model.beta[i, j] <= upperbound[j]

res.__model__.beta_constraint_rule = Constraint(res.__model__.I,
                                                res.__model__.J,
                                                rule=constraint_rule,
                                                doc='beta constraint')

A simple CQR example:

# import packages
from pystoned import CQER
from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
from pystoned import dataset as dataset
from pyomo.environ import Constraint

# import the GHG example data
data = dataset.load_GHG_abatement_cost(x_select=['HRSN', 'CPNK', 'GHG'], y_select=['VALK'])

# calculate the quantile model
res = CQER.CQR(y=data.y, x=data.x, tau=0.5, z=None, cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
# add an extral constraint on beta
def constraint_rule(model, i, j):
    upperbound = [1.0, 1.0, 1.0]
    return model.beta[i, j] <= upperbound[j]

res.__model__.beta_constraint_rule = Constraint(res.__model__.I,
                                                res.__model__.J,
                                                rule=constraint_rule,
                                                doc='beta constraint')
res.optimize(OPT_LOCAL)

# display estimated beta
res.display_beta()

Miscellaneous

gams2python

pyomo provides a good coding environment that can help us smoothly transfer from GAMS to Python. Thus, we prepare a short tutorial to help GAMSers to understand how to rewrite the CNLS models, even other complicated models on Python.

A short comparison:

### We first present the GAMS code,
#sets 
#     i        "DMU's"  /1*10/
#     j        'outputs' /Energy, Length, Customers/;
    
### then the Python code is added to compare.
from pyomo.environ import *
model = ConcreteModel(name = "CNLS")
model.i = Set(initialize=['i1', 'i2', 'i3', 'i4', 'i5', 'i6','i7', 'i8', 'i9','i10'], doc='DMUS', ordered=True)
model.j = Set(initialize=['Energy', 'Length', 'Customers'], doc='outputs')

# alias(i,h); 

model.h = SetOf(model.i)

#Table data(i,j)
#$Include energy.txt;

dtab = {
        ('i1',  'Energy')   : 75,
        ('i1',  'Length')   : 878,
        ('i1',  'Customers'): 4933,
        ('i2',  'Energy')   : 62,
        ('i2',  'Length')   : 964,
        ('i2',  'Customers'): 6149,
        ('i3',  'Energy')   : 78,
        ('i3',  'Length')   : 676,
        ('i3',  'Customers'): 6098,
        ('i4',  'Energy')   : 683,
        ('i4',  'Length')   : 12522,
        ('i4',  'Customers'): 55226,
        ('i5',  'Energy')   : 27,
        ('i5',  'Length')   : 697,
        ('i5',  'Customers'): 1670,
        ('i6',  'Energy')   : 295,
        ('i6',  'Length')   : 953,
        ('i6',  'Customers'): 22949,
        ('i7',  'Energy')   : 44,
        ('i7',  'Length')   : 917,
        ('i7',  'Customers'): 3599,
        ('i8',  'Energy')   : 171,
        ('i8',  'Length')   : 1580,
        ('i8',  'Customers'): 11081,
        ('i9',  'Energy')   : 98,
        ('i9',  'Length')   : 116,
        ('i9',  'Customers'): 377,
        ('i10', 'Energy')   : 203,
        ('i10', 'Length')   : 740,
        ('i10', 'Customers'): 10134,
        }
model.d = Param(model.i, model.j, initialize=dtab, doc='output data')

#PARAMETERS
#c(i)         "Total cost of firm i"
#y(i,j)       "Output j of firm i";
model.c = Param(model.i, initialize={'i1': 1612,
                                     'i2': 1659,
                                     'i3': 1708,
                                     'i4': 18918,
                                     'i5': 1167,
                                     'i6': 3395,
                                     'i7': 1333,
                                     'i8': 3518,
                                     'i9': 1415,
                                     'i10':2469,
                                    }, 
                      doc='Cost data')

def y_init(model, i, j):
    return  model.d[i, j]
model.y = Param(model.i, model.j, initialize=y_init, doc='output data')

#VARIABLES
#E(i)            "Composite error term (v + u)"
#SSE             "Sum of squares of residuals";
#POSITIVE VARIABLES
#b(i,j)    "Beta-coefficients (positivity ensures monotonicity)"
#Chat(i)  ;

model.b = Var(model.i, model.j, bounds=(0.0,None), doc='beta-coeff')
model.e = Var(model.i, doc='res')
model.f = Var(model.i, bounds=(0.0,None), doc='frontier')

#Equations
#QSSE                  objective function = sum of squares of residuals
#QREGRESSION(i)        log-transformed regression equation
#Qlog(i)               supporting hyperplanes of the nonparametric cost function
#QCONC(i,h)            concavity constraint (Afriat inequalities);

#QSSE..                SSE=e=sum(i,E(i)*E(i)) ;
#QREGRESSION(i)..      log(C(i)) =e= log(Chat(i) + 1) + E(i);
#Qlog(i)..             Chat(i) =e= sum(j, b(i,j)*Y(i,j)) - 1;
#QCONC(i,h)..          sum(j, b(i,j)*Y(i,j)) =g= sum(j, b(h,j)*Y(i,j));

def objective_rule(model):
    return sum(model.e[i]*model.e[i] for i in model.i)
model.objective = Objective(rule=objective_rule, sense=minimize, doc='objective function')

def qreg_rule(model, i):
    return log(model.c[i]) == log(model.f[i] + 1) + model.e[i]
model.qreg = Constraint(model.i, rule=qreg_rule, doc='log-transformed regression')

def qlog_rule(model, i):
    return model.f[i] == sum(model.b[i, j]*model.y[i, j] for j in model.j) - 1
model.qlog = Constraint(model.i, rule=qlog_rule, doc='cost function')

def qconvex_rule(model, i, h):
    return sum(model.b[i,j]*model.y[i,j] for j in model.j) >= sum(model.b[h,j]*model.y[i,j] for j in model.j)
model.qconvex = Constraint(model.i, model.h, rule=qconvex_rule, doc='convexity constraint')

# Execute the model
#MODEL StoNED /all/;
#SOLVE StoNED using NLP Minimizing SSE;
    
from pyomo.opt import SolverFactory
import pyomo.environ
from os import environ
solver_manager = SolverManagerFactory('neos')
environ['NEOS_EMAIL'] = 'email@address'
results = solver_manager.solve(model, opt='minos')

#display E.l, b.l;

model.e.display()
model.b.display()      

CNLS_ConcreteModel

We also prepare a concrete model that does not need to call any of our developed functions in the pyStoNED package. In this concrete model, one can define the parameter specification by reading data from an Excel file.

A concreteModel for CNLS estimation:

# Import packages
import pd as pd
import numpy as np
import sys
from pyomo.environ import *

# Sets
i = np.array(['i{}'.format(i) for i in range(1, 90)])
model.i = Set(initialize=i, doc='DMUs', ordered=True )
model.j = Set(initialize=['j1', 'j2', 'j3'], doc='inputs and outputs')

# alias
model.h = SetOf(model.i) 

# Parameters 
# output (y1,y2,y3), avaiable at https://raw.githubusercontent.com/ds2010/pyStoNED/master/pystoned/data/electricityFirms.csv.
df1 = pd.read_excel('y.xlsx', header=0, index_col=0)
Dict1 = dict()
for i in df1.index:
    for j in df1.columns:
        Dict1[i, j] = float(df1[j][i])
model.y = Param(model.i, model.j, initialize = Dict1) 

# input (cost)
df2 = pd.read_excel('x.xlsx', header=0, index_col=0)
Dict2 = dict()
for i in df2.index:
    Dict2[i] = float(df2['TOTEX'][i])
model.c = Param(model.i, initialize = Dict2) 

#Variables
model.b = Var(model.i, model.j, bounds=(0.0,None), doc='beta-coeff')
model.e = Var(model.i, doc='res')
model.f = Var(model.i, bounds=(0.0,None), doc='frontier')

# Constraints and objective f.
def qreg_rule(model, i):
    return log(model.c[i]) == log(model.f[i] + 1) + model.e[i]
model.qreg = Constraint(model.i, rule=qreg_rule, doc='log-transformed regression')

def qlog_rule(model, i):
    return model.f[i] == sum(model.b[i, j]*model.y[i, j] for j in model.j) - 1
model.qlog = Constraint(model.i, rule=qlog_rule, doc='cost function')

def qconcav_rule(model, i, h):
    return sum(model.b[i,j]*model.y[i,j] for j in model.j) >= sum(model.b[h,j]*model.y[i,j] for j in model.j)
model.qconcav = Constraint(model.i, model.h, rule=qconcav_rule, doc='concavity constraint')

def objective_rule(model):
    return sum(model.e[i]*model.e[i] for i in model.i)
model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')

#Solve the model
from os import environ
solver_manager = SolverManagerFactory('neos')
environ['NEOS_EMAIL'] = 'email@address'
results = solver_manager.solve(model, opt='knitro', tee=True)

#retrive the beta
ind = list(model.b)
val = list(model.b[:,:].value)
df_b = [ i + tuple([j]) for i, j in zip(ind, val)]
df_b = np.asarray(df_b)
# display beta coefficients
df_b

df_b = pd.DataFrame(df_b,columns = ['Name', 'Key', 'Value'])
df_b = df_b.pivot(index='Name', columns='Key', values='Value')
# retrive the residuals
ind = list(model.e)
val = list(model.e[:].value)
df_e = np.asarray(val)
#display residuals
df_e

DEA_ConcreteModel

We also prepare a concrete model that can be used to calculate the input oriented VRS model.

A ConcreteModel for DEA estimation:

# import PYOMO package
from pyomo.environ import *
# creat a concrete model
model = ConcreteModel()

# Sets
model.i = Set(initialize=['i1', 'i2','i3', 'i4','i5', 'i6','i7', 'i8','i9', 'i10'], doc='DMUs', ordered=True)
model.j = Set(initialize=['OPEX', 'CAPEX'])   # inputs
model.k = Set(initialize=['Energy', 'Length', 'Customers'])   #outputs

# Alias
model.io = SetOf(model.i) 

# Parameters (define and import the input-output data)
inputdata = {
                ('OPEX', 'i1') :   681,
                ('OPEX', 'i2') :   559,
                ('OPEX', 'i3') :   836,
                ('OPEX', 'i4') :   7559,
                ('OPEX', 'i5') :   424,
                ('OPEX', 'i6') :   1483,
                ('OPEX', 'i7') :   658,
                ('OPEX', 'i8') :   1433,
                ('OPEX', 'i9') :   850,
                ('OPEX', 'i10') :  1155,
                ('CAPEX', 'i1') :  729,
                ('CAPEX', 'i2') :  673,
                ('CAPEX', 'i3') :  851,
                ('CAPEX', 'i4') :  8384,
                ('CAPEX', 'i5') :  562,
                ('CAPEX', 'i6') :  1587,
                ('CAPEX', 'i7') :  570,
                ('CAPEX', 'i8') :  1311,
                ('CAPEX', 'i9') :  564,
                ('CAPEX', 'i10'):  1108,
                }
model.x = Param(model.j, model.i, initialize=inputdata)

outputdata = {
                ('Energy', 'i1'):   75,
                ('Energy', 'i2'):   62,
                ('Energy', 'i3'):   78,
                ('Energy', 'i4'):   683,
                ('Energy', 'i5'):   27,
                ('Energy', 'i6'):   295,
                ('Energy', 'i7'):   44,
                ('Energy', 'i8'):   171,
                ('Energy', 'i9'):   98,
                ('Energy', 'i10'):  203,
                ('Length', 'i1'):   878,
                ('Length', 'i2'):   964,
                ('Length', 'i3'):   676,
                ('Length', 'i4'):   12522,
                ('Length', 'i5'):   697,
                ('Length', 'i6'):   953,
                ('Length', 'i7'):   917,
                ('Length', 'i8'):   1580,
                ('Length', 'i9'):   116,
                ('Length', 'i10'):  740,
                ('Customers', 'i1'):4933,
                ('Customers', 'i2'):6149,
                ('Customers', 'i3'):6098,
                ('Customers', 'i4'):55226,
                ('Customers', 'i5'):1670,
                ('Customers', 'i6'):22949,
                ('Customers', 'i7'):3599,
                ('Customers', 'i8'):11081,
                ('Customers', 'i9'):377,
                ('Customers', 'i10'):10134,
                }
model.y = Param(model.k, model.i, initialize=outputdata)  

# Variables
model.lamda = Var(model.io, model.i, bounds=(0.0, None), doc='envelopment efficiency') 
model.theta = Var(model.io, doc='intensity variable') 

# objective
def objective_rule(model):
    return sum(model.theta[io] for io in model.io)
model.objective = Objective(rule=objective_rule, sense=minimize, doc='Define objective function')       
# Constraints
def output_rule(model, io, k):
    return sum(model.lamda[io, i] * model.y[k, i] for i in model.i) >= model.y[k, io]
model.output = Constraint(model.io, model.k, rule=output_rule, doc='output constraints')
def input_rule(model, io, j):
    return (model.theta[io] * model.x[j, io]) >= sum(
             model.lamda[io, i] * model.x[j, i] for i in model.i)
model.input = Constraint(model.io, model.j, rule=input_rule, doc='input constraints')
def vrs_rule(model, io):
    return sum(model.lamda[io, i] for i in model.i) == 1
model.vrs = Constraint(model.io, rule=vrs_rule, doc='VRS constraints')

# calculate the DEA model 
from pyomo.opt import SolverFactory
from os import environ
solver_manager = SolverManagerFactory('neos')
environ['NEOS_EMAIL'] = 'email@address'
results = solver_manager.solve(model, opt='cplex')

# display the estimates
# efficiency
model.theta.display()
# intensity
model.lamda.display()