from pyomo.environ import *
import pandas as pd

comp_centers = 4
fac_centers = 2
dist_centers = 3

D = [1200, 1000, 900]
QF = [2000, 2000]
QC = [2000, 2000, 2000, 2000]
a = [[30, 100],[70, 90],[70, 90],[110, 40]]
b = [[40, 150, 160],[210, 200, 180]]
s = [800, 830, 880]
c = [80, 70, 80, 50]
pC = [50, 30, 30, 20]
e = [140, 120]
n = [0.65, 0.8, 0.8, 0.9]
m1 = [0.65, 0.9]
r = [0.65, 0.7, 0.75]
lbC = [[160, 230],[170, 190],[180, 200],[180, 110]]
ubC = [[340, 530],[350, 490],[360, 500],[360, 410]]
lbF = [[520, 630, 640],[630, 620, 600]]
ubF = [[800, 830, 880],[800, 830, 880]]
f_fc = 1.0
f_df = 1.0

m = ConcreteModel()

# Variables
m.btp_comp = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.btp_fac = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.btp_dist = Var(RangeSet(dist_centers),domain=NonNegativeReals)
m.btl_comp = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.btl_fac = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.btl_dist = Var(RangeSet(dist_centers),domain=NonNegativeReals)
m.g = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.h = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.x = Var(RangeSet(comp_centers),RangeSet(fac_centers),domain=NonNegativeIntegers)
m.y = Var(RangeSet(fac_centers),RangeSet(dist_centers),domain=NonNegativeIntegers)
m.TP_fc = Var(RangeSet(comp_centers),RangeSet(fac_centers),domain=NonNegativeReals)
m.TP_df = Var(RangeSet(fac_centers),RangeSet(dist_centers),domain=NonNegativeReals)

# Constraints
# Capacity Constraints
m.cap_cons_comp = ConstraintList()
for i in range(comp_centers):
    m.cap_cons_comp.add(expr=m.g[i+1]<=QC[i])

m.cap_cons_fac = ConstraintList()
for j in range(fac_centers):
    m.cap_cons_fac.add(expr=m.h[j+1]<=QF[j])
    
# Balance constraints
m.bal_comp = ConstraintList()
for i in range(comp_centers):
    m.bal_comp.add(expr=sum(m.x[i+1,j+1] for j in range(fac_centers))==m.g[i+1])
    
m.bal_fac = ConstraintList()
for j in range(fac_centers):
    m.bal_comp.add(expr=sum(m.y[j+1,k+1] for k in range(dist_centers))==m.h[j+1])

# Intermediate products and final products balance constraints
m.int_final_bal = ConstraintList()
for j in range(fac_centers):
    m.int_final_bal.add(expr=m.x[1,j+1]+m.x[2,j+1] == sum(m.y[j+1,k+1] for k in range(dist_centers))) # Component 1
for j in range(fac_centers):
    m.int_final_bal.add(expr=m.x[3,j+1]+m.x[4,j+1] == sum(m.y[j+1,k+1] for k in range(dist_centers))) # Component 2
    
# Demand satisfaction
m.dem_cons = ConstraintList()
for k in range(dist_centers):
    m.dem_cons.add(expr=sum(m.y[j+1,k+1] for j in range(fac_centers))==D[k])
    
# Transfer prices bounds
m.tp_bnds_comp = ConstraintList()
for i in range(comp_centers):
    for j in range(fac_centers):
        m.tp_bnds_comp.add(expr=lbC[i][j]*m.x[i+1,j+1] <= m.TP_fc[i+1,j+1])
    
for i in range(comp_centers):
    for j in range(fac_centers):
        m.tp_bnds_comp.add(expr=m.TP_fc[i+1,j+1] <= ubC[i][j]*m.x[i+1,j+1])

m.tp_bnds_fac = ConstraintList()
for j in range(fac_centers):
    for k in range(dist_centers):
        m.tp_bnds_fac.add(expr=lbF[j][k]*m.y[j+1,k+1] <= m.TP_df[j+1,k+1])

for j in range(fac_centers):
    for k in range(dist_centers):
        m.tp_bnds_fac.add(expr=m.TP_df[j+1,k+1] <= ubF[j][k]*m.y[j+1,k+1])

# Profit before tax computation constraints
m.comp_btp_cons = ConstraintList() # component production centers profit before tax
for i in range(comp_centers):
    m.comp_btp_cons.add(expr=m.btp_comp[i+1]-m.btl_comp[i+1] == sum((m.TP_fc[i+1,j+1]-f_fc*a[i][j]*m.x[i+1,j+1])for j in range(fac_centers))-m.g[i+1]*(pC[i]+c[i]))
                        
m.fac_btp_cons = ConstraintList() # Final product production facilities
for j in range(fac_centers):
    m.fac_btp_cons.add(expr=m.btp_fac[j+1]-m.btl_fac[j+1] == sum((m.TP_df[j+1,k+1]-f_df*b[j][k]*m.y[j+1,k+1]) for k in range(dist_centers)) - e[j]*m.h[j+1] - sum((m.TP_fc[i+1,j+1]+(1-f_fc)*a[i][j]*m.x[i+1,j+1]) for i in range(comp_centers)))
    
m.dist_btp_cons = ConstraintList()
for k in range(dist_centers):
    m.dist_btp_cons.add(expr= m.btp_dist[k+1]-m.btl_dist[k+1] == sum((s[k]*m.y[j+1,k+1]-m.TP_df[j+1,k+1]-(1-f_df)*b[j][k]*m.y[j+1,k+1])for j in range(fac_centers)))
    
# Objective
#m.obj = Objective(expr = sum((1-n[i])*m.btp_comp[i+1]-m.btl_comp[i+1] for i in range(comp_centers))+sum((1-m1[j])*m.btp_fac[j+1]-m.btl_fac[j+1] for j in range(fac_centers))+sum((1-r[k])*m.btp_dist[k+1]-m.btl_dist[k+1] for k in range(dist_centers)),sense=maximize)

m.obj = Objective(expr = sum((n[i])*m.btp_comp[i+1]-m.btl_comp[i+1] for i in range(comp_centers))+sum((m1[j])*m.btp_fac[j+1]-m.btl_fac[j+1] for j in range(fac_centers))+sum((r[k])*m.btp_dist[k+1]-m.btl_dist[k+1] for k in range(dist_centers)),sense=maximize)

# Solving the model
opt = SolverFactory('cplex')
results = opt.solve(m,tee=True)
print("Solver status:", results.solver.status)
print("Solver termination condition:", results.solver.termination_condition)


cols = ['Variable','Variable Name','Number','TP_variable','TP_fixed_mean','TP_fixed_min','TP_fixed_max']
sol = pd.DataFrame(data=None,columns=cols,index=range(35))

ind = 0

for i in range(comp_centers):
    sol.at[ind,'Variable'] = 'Comp Prod'
    sol.at[ind,'Variable Name'] = 'g'
    sol.at[ind,'Number'] = i+1
    sol.at[ind,'TP_variable'] = m.g[i+1].value
    ind += 1
    
for i in range(fac_centers):
    sol.at[ind,'Variable'] = 'Final Prod'
    sol.at[ind,'Variable Name'] = 'h'
    sol.at[ind,'Number'] = i+1
    sol.at[ind,'TP_variable'] = m.h[i+1].value
    ind += 1

for i in range(comp_centers):
    for j in range(fac_centers):
        sol.at[ind,'Variable'] = 'Comp ship'
        sol.at[ind,'Number'] = (i+1,j+1)
        sol.at[ind,'Variable Name'] = 'x'
        sol.at[ind,'TP_variable'] = m.x[i+1,j+1].value
        ind += 1

for i in range(fac_centers):
    for j in range(dist_centers):
        sol.at[ind,'Variable'] = 'Prod ship'
        sol.at[ind,'Number'] = (i+1,j+1)
        sol.at[ind,'Variable Name'] = 'y'
        sol.at[ind,'TP_variable'] = m.y[i+1,j+1].value
        ind += 1

for i in range(comp_centers):
    for j in range(fac_centers):
        sol.at[ind,'Variable'] = 'TP c,f'
        sol.at[ind,'Number'] = (i+1,j+1)
        sol.at[ind,'Variable Name'] = 'tp'
        if m.TP_fc[i+1,j+1].value == 0 and m.x[i+1,j+1].value == 0:
            sol.at[ind,'TP_variable'] = 0
        else:
            sol.at[ind,'TP_variable'] = m.TP_fc[i+1,j+1].value/m.x[i+1,j+1].value
        ind += 1
    
for i in range(fac_centers):
    for j in range(dist_centers):
        sol.at[ind,'Variable'] = 'TP f,d'
        sol.at[ind,'Number'] = (i+1,j+1)
        sol.at[ind,'Variable Name'] = 'tp'
        if m.TP_df[i+1,j+1].value == 0 and m.y[i+1,j+1].value == 0:
            sol.at[ind,'TP_variable'] = 0
        else:
            sol.at[ind,'TP_variable'] = m.TP_df[i+1,j+1].value/m.y[i+1,j+1].value
        ind += 1
        
sol.at[ind,'Variable'] = 'Profit'
sol.at[ind,'TP_variable'] = m.obj() 
ind += 1
sol.at[ind,'Variable'] = 'Comp Man Profit'
sol.at[ind,'TP_variable'] = sum((n[i])*m.btp_comp[i+1].value-m.btl_comp[i+1].value for i in range(comp_centers))
ind += 1
sol.at[ind,'Variable'] = 'Prod Man Profit'
sol.at[ind,'TP_variable'] = sum((m1[j])*m.btp_fac[j+1].value-m.btl_fac[j+1].value for j in range(fac_centers))
ind += 1
sol.at[ind,'Variable'] = 'Dist Profit'
sol.at[ind,'TP_variable'] = sum((r[k])*m.btp_dist[k+1].value-m.btl_dist[k+1].value for k in range(dist_centers))
##########################################################
##########################################################
##########################################################

m = ConcreteModel()

# Variables
m.btp_comp = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.btp_fac = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.btp_dist = Var(RangeSet(dist_centers),domain=NonNegativeReals)
m.btl_comp = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.btl_fac = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.btl_dist = Var(RangeSet(dist_centers),domain=NonNegativeReals)
m.g = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.h = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.x = Var(RangeSet(comp_centers),RangeSet(fac_centers),domain=NonNegativeIntegers)
m.y = Var(RangeSet(fac_centers),RangeSet(dist_centers),domain=NonNegativeIntegers)
m.TP_fc = Var(RangeSet(comp_centers),RangeSet(fac_centers),domain=NonNegativeReals)
m.TP_df = Var(RangeSet(fac_centers),RangeSet(dist_centers),domain=NonNegativeReals)

# Constraints
# Capacity Constraints
m.cap_cons_comp = ConstraintList()
for i in range(comp_centers):
    m.cap_cons_comp.add(expr=m.g[i+1]<=QC[i])

m.cap_cons_fac = ConstraintList()
for j in range(fac_centers):
    m.cap_cons_fac.add(expr=m.h[j+1]<=QF[j])
    
# Balance constraints
m.bal_comp = ConstraintList()
for i in range(comp_centers):
    m.bal_comp.add(expr=sum(m.x[i+1,j+1] for j in range(fac_centers))==m.g[i+1])
    
m.bal_fac = ConstraintList()
for j in range(fac_centers):
    m.bal_comp.add(expr=sum(m.y[j+1,k+1] for k in range(dist_centers))==m.h[j+1])

# Intermediate products and final products balance constraints
m.int_final_bal = ConstraintList()
for j in range(fac_centers):
    m.int_final_bal.add(expr=m.x[1,j+1]+m.x[2,j+1] == sum(m.y[j+1,k+1] for k in range(dist_centers))) # Component 1
for j in range(fac_centers):
    m.int_final_bal.add(expr=m.x[3,j+1]+m.x[4,j+1] == sum(m.y[j+1,k+1] for k in range(dist_centers))) # Component 2
    
# Demand satisfaction
m.dem_cons = ConstraintList()
for k in range(dist_centers):
    m.dem_cons.add(expr=sum(m.y[j+1,k+1] for j in range(fac_centers))==D[k])
    
# Transfer prices bounds
m.tp_bnds_comp = ConstraintList()
for i in range(comp_centers):
    for j in range(fac_centers):
        m.tp_bnds_comp.add(expr=(lbC[i][j]+ubC[i][j])/2*m.x[i+1,j+1] <= m.TP_fc[i+1,j+1])
    
#for i in range(comp_centers):
#    for j in range(fac_centers):
#        m.tp_bnds_comp.add(expr=m.TP_fc[i+1,j+1] <= ubC[i][j]*m.x[i+1,j+1])

m.tp_bnds_fac = ConstraintList()
for j in range(fac_centers):
    for k in range(dist_centers):
        m.tp_bnds_fac.add(expr=(lbF[j][k]+ubF[j][k])/2*m.y[j+1,k+1] == m.TP_df[j+1,k+1])

#for j in range(fac_centers):
#    for k in range(dist_centers):
#        m.tp_bnds_fac.add(expr=m.TP_df[j+1,k+1] <= ubF[j][k]*m.y[j+1,k+1])

# Profit before tax computation constraints
m.comp_btp_cons = ConstraintList() # component production centers profit before tax
for i in range(comp_centers):
    m.comp_btp_cons.add(expr=m.btp_comp[i+1]-m.btl_comp[i+1] == sum((m.TP_fc[i+1,j+1]-f_fc*a[i][j]*m.x[i+1,j+1])for j in range(fac_centers))-m.g[i+1]*(pC[i]+c[i]))
                        
m.fac_btp_cons = ConstraintList() # Final product production facilities
for j in range(fac_centers):
    m.fac_btp_cons.add(expr=m.btp_fac[j+1]-m.btl_fac[j+1] == sum((m.TP_df[j+1,k+1]-f_df*b[j][k]*m.y[j+1,k+1]) for k in range(dist_centers)) - e[j]*m.h[j+1] - sum((m.TP_fc[i+1,j+1]+(1-f_fc)*a[i][j]*m.x[i+1,j+1]) for i in range(comp_centers)))
    
m.dist_btp_cons = ConstraintList()
for k in range(dist_centers):
    m.dist_btp_cons.add(expr= m.btp_dist[k+1]-m.btl_dist[k+1] == sum((s[k]*m.y[j+1,k+1]-m.TP_df[j+1,k+1]-(1-f_df)*b[j][k]*m.y[j+1,k+1])for j in range(fac_centers)))
    
# Objective
#m.obj = Objective(expr = sum((1-n[i])*m.btp_comp[i+1]-m.btl_comp[i+1] for i in range(comp_centers))+sum((1-m1[j])*m.btp_fac[j+1]-m.btl_fac[j+1] for j in range(fac_centers))+sum((1-r[k])*m.btp_dist[k+1]-m.btl_dist[k+1] for k in range(dist_centers)),sense=maximize)

m.obj = Objective(expr = sum((n[i])*m.btp_comp[i+1]-m.btl_comp[i+1] for i in range(comp_centers))+sum((m1[j])*m.btp_fac[j+1]-m.btl_fac[j+1] for j in range(fac_centers))+sum((r[k])*m.btp_dist[k+1]-m.btl_dist[k+1] for k in range(dist_centers)),sense=maximize)

# Solving the model
opt = SolverFactory('cplex')
results = opt.solve(m,tee=True)
print("Solver status:", results.solver.status)
print("Solver termination condition:", results.solver.termination_condition)


ind = 0

for i in range(comp_centers):
    #sol.at[ind,'Variable'] = 'Comp Prod'
    #sol.at[ind,'Variable Name'] = 'g'
    #sol.at[ind,'Number'] = i+1
    sol.at[ind,'TP_fixed_mean'] = m.g[i+1].value
    ind += 1
    
for i in range(fac_centers):
    #sol.at[ind,'Variable'] = 'Final Prod'
    #sol.at[ind,'Variable Name'] = 'h'
    #sol.at[ind,'Number'] = i+1
    sol.at[ind,'TP_fixed_mean'] = m.h[i+1].value
    ind += 1

for i in range(comp_centers):
    for j in range(fac_centers):
        #sol.at[ind,'Variable'] = 'Comp ship'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'x'
        sol.at[ind,'TP_fixed_mean'] = m.x[i+1,j+1].value
        ind += 1

for i in range(fac_centers):
    for j in range(dist_centers):
        #sol.at[ind,'Variable'] = 'Prod ship'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'y'
        sol.at[ind,'TP_fixed_mean'] = m.y[i+1,j+1].value
        ind += 1

for i in range(comp_centers):
    for j in range(fac_centers):
        #sol.at[ind,'Variable'] = 'TP c,f'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'tp'
        if m.TP_fc[i+1,j+1].value == 0 and m.x[i+1,j+1].value == 0:
            sol.at[ind,'TP_fixed_mean'] = 0
        else:
            sol.at[ind,'TP_fixed_mean'] = m.TP_fc[i+1,j+1].value/m.x[i+1,j+1].value
        ind += 1
    
for i in range(fac_centers):
    for j in range(dist_centers):
        #sol.at[ind,'Variable'] = 'TP f,d'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'tp'
        if m.TP_df[i+1,j+1].value == 0 and m.y[i+1,j+1].value == 0:
            sol.at[ind,'TP_fixed_mean'] = 0
        else:
            sol.at[ind,'TP_fixed_mean'] = m.TP_df[i+1,j+1].value/m.y[i+1,j+1].value
        ind += 1
        
#sol.at[ind,'Variable'] = 'Profit'
sol.at[ind,'TP_fixed_mean'] = m.obj() 
ind += 1
#sol.at[ind,'Variable'] = 'Comp Man Profit'
sol.at[ind,'TP_fixed_mean'] = sum((n[i])*m.btp_comp[i+1].value-m.btl_comp[i+1].value for i in range(comp_centers))
ind += 1
#sol.at[ind,'Variable'] = 'Prod Man Profit'
sol.at[ind,'TP_fixed_mean'] = sum((m1[j])*m.btp_fac[j+1].value-m.btl_fac[j+1].value for j in range(fac_centers))
ind += 1
#sol.at[ind,'Variable'] = 'Dist Profit'
sol.at[ind,'TP_fixed_mean'] = sum((r[k])*m.btp_dist[k+1].value-m.btl_dist[k+1].value for k in range(dist_centers))

########################################################
########################################################
########################################################
########################################################

m = ConcreteModel()

# Variables
m.btp_comp = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.btp_fac = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.btp_dist = Var(RangeSet(dist_centers),domain=NonNegativeReals)
m.btl_comp = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.btl_fac = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.btl_dist = Var(RangeSet(dist_centers),domain=NonNegativeReals)
m.g = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.h = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.x = Var(RangeSet(comp_centers),RangeSet(fac_centers),domain=NonNegativeIntegers)
m.y = Var(RangeSet(fac_centers),RangeSet(dist_centers),domain=NonNegativeIntegers)
m.TP_fc = Var(RangeSet(comp_centers),RangeSet(fac_centers),domain=NonNegativeReals)
m.TP_df = Var(RangeSet(fac_centers),RangeSet(dist_centers),domain=NonNegativeReals)

# Constraints
# Capacity Constraints
m.cap_cons_comp = ConstraintList()
for i in range(comp_centers):
    m.cap_cons_comp.add(expr=m.g[i+1]<=QC[i])

m.cap_cons_fac = ConstraintList()
for j in range(fac_centers):
    m.cap_cons_fac.add(expr=m.h[j+1]<=QF[j])
    
# Balance constraints
m.bal_comp = ConstraintList()
for i in range(comp_centers):
    m.bal_comp.add(expr=sum(m.x[i+1,j+1] for j in range(fac_centers))==m.g[i+1])
    
m.bal_fac = ConstraintList()
for j in range(fac_centers):
    m.bal_comp.add(expr=sum(m.y[j+1,k+1] for k in range(dist_centers))==m.h[j+1])

# Intermediate products and final products balance constraints
m.int_final_bal = ConstraintList()
for j in range(fac_centers):
    m.int_final_bal.add(expr=m.x[1,j+1]+m.x[2,j+1] == sum(m.y[j+1,k+1] for k in range(dist_centers))) # Component 1
for j in range(fac_centers):
    m.int_final_bal.add(expr=m.x[3,j+1]+m.x[4,j+1] == sum(m.y[j+1,k+1] for k in range(dist_centers))) # Component 2
    
# Demand satisfaction
m.dem_cons = ConstraintList()
for k in range(dist_centers):
    m.dem_cons.add(expr=sum(m.y[j+1,k+1] for j in range(fac_centers))==D[k])
    
# Transfer prices bounds
m.tp_bnds_comp = ConstraintList()
for i in range(comp_centers):
    for j in range(fac_centers):
        m.tp_bnds_comp.add(expr=(lbC[i][j])*m.x[i+1,j+1] <= m.TP_fc[i+1,j+1])
    
#for i in range(comp_centers):
#    for j in range(fac_centers):
#        m.tp_bnds_comp.add(expr=m.TP_fc[i+1,j+1] <= ubC[i][j]*m.x[i+1,j+1])

m.tp_bnds_fac = ConstraintList()
for j in range(fac_centers):
    for k in range(dist_centers):
        m.tp_bnds_fac.add(expr=(lbF[j][k])*m.y[j+1,k+1] == m.TP_df[j+1,k+1])

#for j in range(fac_centers):
#    for k in range(dist_centers):
#        m.tp_bnds_fac.add(expr=m.TP_df[j+1,k+1] <= ubF[j][k]*m.y[j+1,k+1])

# Profit before tax computation constraints
m.comp_btp_cons = ConstraintList() # component production centers profit before tax
for i in range(comp_centers):
    m.comp_btp_cons.add(expr=m.btp_comp[i+1]-m.btl_comp[i+1] == sum((m.TP_fc[i+1,j+1]-f_fc*a[i][j]*m.x[i+1,j+1])for j in range(fac_centers))-m.g[i+1]*(pC[i]+c[i]))
                        
m.fac_btp_cons = ConstraintList() # Final product production facilities
for j in range(fac_centers):
    m.fac_btp_cons.add(expr=m.btp_fac[j+1]-m.btl_fac[j+1] == sum((m.TP_df[j+1,k+1]-f_df*b[j][k]*m.y[j+1,k+1]) for k in range(dist_centers)) - e[j]*m.h[j+1] - sum((m.TP_fc[i+1,j+1]+(1-f_fc)*a[i][j]*m.x[i+1,j+1]) for i in range(comp_centers)))
    
m.dist_btp_cons = ConstraintList()
for k in range(dist_centers):
    m.dist_btp_cons.add(expr= m.btp_dist[k+1]-m.btl_dist[k+1] == sum((s[k]*m.y[j+1,k+1]-m.TP_df[j+1,k+1]-(1-f_df)*b[j][k]*m.y[j+1,k+1])for j in range(fac_centers)))
    
# Objective
#m.obj = Objective(expr = sum((1-n[i])*m.btp_comp[i+1]-m.btl_comp[i+1] for i in range(comp_centers))+sum((1-m1[j])*m.btp_fac[j+1]-m.btl_fac[j+1] for j in range(fac_centers))+sum((1-r[k])*m.btp_dist[k+1]-m.btl_dist[k+1] for k in range(dist_centers)),sense=maximize)

m.obj = Objective(expr = sum((n[i])*m.btp_comp[i+1]-m.btl_comp[i+1] for i in range(comp_centers))+sum((m1[j])*m.btp_fac[j+1]-m.btl_fac[j+1] for j in range(fac_centers))+sum((r[k])*m.btp_dist[k+1]-m.btl_dist[k+1] for k in range(dist_centers)),sense=maximize)

# Solving the model
opt = SolverFactory('cplex')
results = opt.solve(m,tee=True)
print("Solver status:", results.solver.status)
print("Solver termination condition:", results.solver.termination_condition)


ind = 0

for i in range(comp_centers):
    #sol.at[ind,'Variable'] = 'Comp Prod'
    #sol.at[ind,'Variable Name'] = 'g'
    #sol.at[ind,'Number'] = i+1
    sol.at[ind,'TP_fixed_min'] = m.g[i+1].value
    ind += 1
    
for i in range(fac_centers):
    #sol.at[ind,'Variable'] = 'Final Prod'
    #sol.at[ind,'Variable Name'] = 'h'
    #sol.at[ind,'Number'] = i+1
    sol.at[ind,'TP_fixed_min'] = m.h[i+1].value
    ind += 1

for i in range(comp_centers):
    for j in range(fac_centers):
        #sol.at[ind,'Variable'] = 'Comp ship'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'x'
        sol.at[ind,'TP_fixed_min'] = m.x[i+1,j+1].value
        ind += 1

for i in range(fac_centers):
    for j in range(dist_centers):
        #sol.at[ind,'Variable'] = 'Prod ship'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'y'
        sol.at[ind,'TP_fixed_min'] = m.y[i+1,j+1].value
        ind += 1

for i in range(comp_centers):
    for j in range(fac_centers):
        #sol.at[ind,'Variable'] = 'TP c,f'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'tp'
        if m.TP_fc[i+1,j+1].value == 0 and m.x[i+1,j+1].value == 0:
            sol.at[ind,'TP_fixed_min'] = 0
        else:
            sol.at[ind,'TP_fixed_min'] = m.TP_fc[i+1,j+1].value/m.x[i+1,j+1].value
        ind += 1
    
for i in range(fac_centers):
    for j in range(dist_centers):
        #sol.at[ind,'Variable'] = 'TP f,d'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'tp'
        if m.TP_df[i+1,j+1].value == 0 and m.y[i+1,j+1].value == 0:
            sol.at[ind,'TP_fixed_min'] = 0
        else:
            sol.at[ind,'TP_fixed_min'] = m.TP_df[i+1,j+1].value/m.y[i+1,j+1].value
        ind += 1
        
#sol.at[ind,'Variable'] = 'Profit'
sol.at[ind,'TP_fixed_min'] = m.obj() 
ind += 1
#sol.at[ind,'Variable'] = 'Comp Man Profit'
sol.at[ind,'TP_fixed_min'] = sum((n[i])*m.btp_comp[i+1].value-m.btl_comp[i+1].value for i in range(comp_centers))
ind += 1
#sol.at[ind,'Variable'] = 'Prod Man Profit'
sol.at[ind,'TP_fixed_min'] = sum((m1[j])*m.btp_fac[j+1].value-m.btl_fac[j+1].value for j in range(fac_centers))
ind += 1
#sol.at[ind,'Variable'] = 'Dist Profit'
sol.at[ind,'TP_fixed_min'] = sum((r[k])*m.btp_dist[k+1].value-m.btl_dist[k+1].value for k in range(dist_centers))


########################################################
########################################################
########################################################
########################################################

m = ConcreteModel()

# Variables
m.btp_comp = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.btp_fac = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.btp_dist = Var(RangeSet(dist_centers),domain=NonNegativeReals)
m.btl_comp = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.btl_fac = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.btl_dist = Var(RangeSet(dist_centers),domain=NonNegativeReals)
m.g = Var(RangeSet(comp_centers),domain=NonNegativeReals)
m.h = Var(RangeSet(fac_centers),domain=NonNegativeReals)
m.x = Var(RangeSet(comp_centers),RangeSet(fac_centers),domain=NonNegativeIntegers)
m.y = Var(RangeSet(fac_centers),RangeSet(dist_centers),domain=NonNegativeIntegers)
m.TP_fc = Var(RangeSet(comp_centers),RangeSet(fac_centers),domain=NonNegativeReals)
m.TP_df = Var(RangeSet(fac_centers),RangeSet(dist_centers),domain=NonNegativeReals)

# Constraints
# Capacity Constraints
m.cap_cons_comp = ConstraintList()
for i in range(comp_centers):
    m.cap_cons_comp.add(expr=m.g[i+1]<=QC[i])

m.cap_cons_fac = ConstraintList()
for j in range(fac_centers):
    m.cap_cons_fac.add(expr=m.h[j+1]<=QF[j])
    
# Balance constraints
m.bal_comp = ConstraintList()
for i in range(comp_centers):
    m.bal_comp.add(expr=sum(m.x[i+1,j+1] for j in range(fac_centers))==m.g[i+1])
    
m.bal_fac = ConstraintList()
for j in range(fac_centers):
    m.bal_comp.add(expr=sum(m.y[j+1,k+1] for k in range(dist_centers))==m.h[j+1])

# Intermediate products and final products balance constraints
m.int_final_bal = ConstraintList()
for j in range(fac_centers):
    m.int_final_bal.add(expr=m.x[1,j+1]+m.x[2,j+1] == sum(m.y[j+1,k+1] for k in range(dist_centers))) # Component 1
for j in range(fac_centers):
    m.int_final_bal.add(expr=m.x[3,j+1]+m.x[4,j+1] == sum(m.y[j+1,k+1] for k in range(dist_centers))) # Component 2
    
# Demand satisfaction
m.dem_cons = ConstraintList()
for k in range(dist_centers):
    m.dem_cons.add(expr=sum(m.y[j+1,k+1] for j in range(fac_centers))==D[k])
    
# Transfer prices bounds
m.tp_bnds_comp = ConstraintList()
for i in range(comp_centers):
    for j in range(fac_centers):
        m.tp_bnds_comp.add(expr=(ubC[i][j])*m.x[i+1,j+1] <= m.TP_fc[i+1,j+1])
    
#for i in range(comp_centers):
#    for j in range(fac_centers):
#        m.tp_bnds_comp.add(expr=m.TP_fc[i+1,j+1] <= ubC[i][j]*m.x[i+1,j+1])

m.tp_bnds_fac = ConstraintList()
for j in range(fac_centers):
    for k in range(dist_centers):
        m.tp_bnds_fac.add(expr=(ubF[j][k])*m.y[j+1,k+1] == m.TP_df[j+1,k+1])

#for j in range(fac_centers):
#    for k in range(dist_centers):
#        m.tp_bnds_fac.add(expr=m.TP_df[j+1,k+1] <= ubF[j][k]*m.y[j+1,k+1])

# Profit before tax computation constraints
m.comp_btp_cons = ConstraintList() # component production centers profit before tax
for i in range(comp_centers):
    m.comp_btp_cons.add(expr=m.btp_comp[i+1]-m.btl_comp[i+1] == sum((m.TP_fc[i+1,j+1]-f_fc*a[i][j]*m.x[i+1,j+1])for j in range(fac_centers))-m.g[i+1]*(pC[i]+c[i]))
                        
m.fac_btp_cons = ConstraintList() # Final product production facilities
for j in range(fac_centers):
    m.fac_btp_cons.add(expr=m.btp_fac[j+1]-m.btl_fac[j+1] == sum((m.TP_df[j+1,k+1]-f_df*b[j][k]*m.y[j+1,k+1]) for k in range(dist_centers)) - e[j]*m.h[j+1] - sum((m.TP_fc[i+1,j+1]+(1-f_fc)*a[i][j]*m.x[i+1,j+1]) for i in range(comp_centers)))
    
m.dist_btp_cons = ConstraintList()
for k in range(dist_centers):
    m.dist_btp_cons.add(expr= m.btp_dist[k+1]-m.btl_dist[k+1] == sum((s[k]*m.y[j+1,k+1]-m.TP_df[j+1,k+1]-(1-f_df)*b[j][k]*m.y[j+1,k+1])for j in range(fac_centers)))
    
# Objective
#m.obj = Objective(expr = sum((1-n[i])*m.btp_comp[i+1]-m.btl_comp[i+1] for i in range(comp_centers))+sum((1-m1[j])*m.btp_fac[j+1]-m.btl_fac[j+1] for j in range(fac_centers))+sum((1-r[k])*m.btp_dist[k+1]-m.btl_dist[k+1] for k in range(dist_centers)),sense=maximize)

m.obj = Objective(expr = sum((n[i])*m.btp_comp[i+1]-m.btl_comp[i+1] for i in range(comp_centers))+sum((m1[j])*m.btp_fac[j+1]-m.btl_fac[j+1] for j in range(fac_centers))+sum((r[k])*m.btp_dist[k+1]-m.btl_dist[k+1] for k in range(dist_centers)),sense=maximize)

# Solving the model
opt = SolverFactory('cplex')
results = opt.solve(m,tee=True)
print("Solver status:", results.solver.status)
print("Solver termination condition:", results.solver.termination_condition)


ind = 0

for i in range(comp_centers):
    #sol.at[ind,'Variable'] = 'Comp Prod'
    #sol.at[ind,'Variable Name'] = 'g'
    #sol.at[ind,'Number'] = i+1
    sol.at[ind,'TP_fixed_max'] = m.g[i+1].value
    ind += 1
    
for i in range(fac_centers):
    #sol.at[ind,'Variable'] = 'Final Prod'
    #sol.at[ind,'Variable Name'] = 'h'
    #sol.at[ind,'Number'] = i+1
    sol.at[ind,'TP_fixed_max'] = m.h[i+1].value
    ind += 1

for i in range(comp_centers):
    for j in range(fac_centers):
        #sol.at[ind,'Variable'] = 'Comp ship'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'x'
        sol.at[ind,'TP_fixed_max'] = m.x[i+1,j+1].value
        ind += 1

for i in range(fac_centers):
    for j in range(dist_centers):
        #sol.at[ind,'Variable'] = 'Prod ship'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'y'
        sol.at[ind,'TP_fixed_max'] = m.y[i+1,j+1].value
        ind += 1

for i in range(comp_centers):
    for j in range(fac_centers):
        #sol.at[ind,'Variable'] = 'TP c,f'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'tp'
        if m.TP_fc[i+1,j+1].value == 0 and m.x[i+1,j+1].value == 0:
            sol.at[ind,'TP_fixed_max'] = 0
        else:
            sol.at[ind,'TP_fixed_max'] = m.TP_fc[i+1,j+1].value/m.x[i+1,j+1].value
        ind += 1
    
for i in range(fac_centers):
    for j in range(dist_centers):
        #sol.at[ind,'Variable'] = 'TP f,d'
        #sol.at[ind,'Number'] = (i+1,j+1)
        #sol.at[ind,'Variable Name'] = 'tp'
        if m.TP_df[i+1,j+1].value == 0 and m.y[i+1,j+1].value == 0:
            sol.at[ind,'TP_fixed_max'] = 0
        else:
            sol.at[ind,'TP_fixed_max'] = m.TP_df[i+1,j+1].value/m.y[i+1,j+1].value
        ind += 1
        
#sol.at[ind,'Variable'] = 'Profit'
sol.at[ind,'TP_fixed_max'] = m.obj() 
ind += 1
#sol.at[ind,'Variable'] = 'Comp Man Profit'
sol.at[ind,'TP_fixed_max'] = sum((n[i])*m.btp_comp[i+1].value-m.btl_comp[i+1].value for i in range(comp_centers))
ind += 1
#sol.at[ind,'Variable'] = 'Prod Man Profit'
sol.at[ind,'TP_fixed_max'] = sum((m1[j])*m.btp_fac[j+1].value-m.btl_fac[j+1].value for j in range(fac_centers))
ind += 1
#sol.at[ind,'Variable'] = 'Dist Profit'
sol.at[ind,'TP_fixed_max'] = sum((r[k])*m.btp_dist[k+1].value-m.btl_dist[k+1].value for k in range(dist_centers))
