-
Notifications
You must be signed in to change notification settings - Fork 0
/
CLCO.py
2361 lines (2145 loc) · 133 KB
/
CLCO.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import os.path
import glob
import pyomo.environ as pyo
from pyomo.opt import TerminationCondition
import csv
from CLCO_Data import CLCO_Data
from sympy import symbols
import matplotlib.pyplot as plt
import math
def utopian(m, lca_midpoint, lca_type):
"""
Finds the utopia and nadir points for the pareto front with NPV and any LCA subcategory
:param m: the model
:param lca_midpoint: the LCA midpoint type
:param lca_type: the LCA type: ALCA, CLCA
:return: two lists containing the utopia and nadir points
"""
# step 1: normalize the objective functions
# solve the model to get optimal objective function values
utopia = []
nadir = []
# only one function can be active at a time
m.Obj2.deactivate()
m.Obj.activate()
model = m
opt = pyo.SolverFactory('gurobi')
print(opt.solve(model)) # keepfiles = True
utopia.append(pyo.value(model.npv[0]))
temp = pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint])
# only one function can be active at a time
m.Obj.deactivate()
m.Obj2.activate()
model = m
opt = pyo.SolverFactory('gurobi')
print(opt.solve(model)) # keepfiles = True
# get the nadir and utopia points
utopia.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint]))
nadir.append(pyo.value(model.npv[0]))
nadir.append(pyo.value(temp))
# add constraints on the anchor points
m.const.add(expr=m.npv[0] >= nadir[0])
m.const.add(expr=m.total_LCA_midpoints[0, lca_type, lca_midpoint] <= nadir[1])
print("scalar numerator", (nadir[0] - utopia[0]), "nadir npv", nadir[0], "utopia npv", utopia[0])
print("scalar denominator", (nadir[1] - utopia[1]), "nadir lca_midpoint", nadir[1], "utopia lca_midpoint",
utopia[1])
if (nadir[1] - utopia[1]) <= .000001:
raise ValueError("No Pareto Front can be generated - no change in denominator")
scalar = abs((nadir[0] - utopia[0]) / (nadir[1] - utopia[1]))
print("scalar", scalar)
# deactivate remaining objective function
m.Obj2.deactivate()
return utopia, nadir
def utopian3D(m, lca_midpoint1, lca_midpoint2, lca_type):
"""
Finds the utopia and nadir points for the pareto front with NPV and any LCA subcategory
:param m: the model
:param lca_midpoint1: the first LCA midpoint type
:param lca_midpoint2: the second LCA midpoint type
:param lca_type: the LCA type: ALCA, CLCA
:return: two lists containing the utopia and nadir points, as well as the x, y, and z points for the pareto front
"""
# step 1: normalize the objective functions
# solve the model to get optimal objective function values
utopia = []
nadir = []
gwp = []
fe = []
npv = []
# only one function can be active at a time
m.Obj3.deactivate()
m.Obj2.deactivate()
m.Obj.activate()
model = m
opt = pyo.SolverFactory('gurobi')
print(opt.solve(model, tee=True)) # keepfiles = True
# print model
print_model(scenario, model, int(pyo.value(model.npv[0])), "TEA")
print_model(scenario, model, int(pyo.value(model.npv[0])), "LCA", midpoint=midpoint)
utopia.append(pyo.value(model.npv[0]))
temp = []
temp1 = [pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint1])]
temp2 = [pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint2])]
npv.append(pyo.value(model.npv[0]))
gwp.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint1]))
fe.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint2]))
# only one function can be active at a time
m.Obj.deactivate()
m.Obj2.activate()
model = m
opt = pyo.SolverFactory('gurobi')
print(opt.solve(model, tee=True)) # keepfiles = True
# print model
print_model(scenario, model, int(pyo.value(model.npv[0])), "TEA")
print_model(scenario, model, int(pyo.value(model.npv[0])), "LCA", midpoint=midpoint)
# get the nadir and utopia points
utopia.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint1]))
temp.append(pyo.value(model.npv[0]))
temp2.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint2]))
npv.append(pyo.value(model.npv[0]))
gwp.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint1]))
fe.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint2]))
# only one function can be active at a time
m.Obj2.deactivate()
m.Obj3.activate()
model = m
opt = pyo.SolverFactory('gurobi')
print(opt.solve(model, tee=True)) # keepfiles = True
# print model
print_model(scenario, model, int(pyo.value(model.npv[0])), "TEA")
print_model(scenario, model, int(pyo.value(model.npv[0])), "LCA", midpoint=midpoint)
# get the nadir and utopia points
utopia.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint2]))
temp.append(pyo.value(model.npv[0]))
temp1.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint2]))
npv.append(pyo.value(model.npv[0]))
gwp.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint1]))
fe.append(pyo.value(model.total_LCA_midpoints[0, lca_type, lca_midpoint2]))
# capture the smallest points and load them into an array
nadir.append(min(temp))
nadir.append(max(temp1))
nadir.append(max(temp2))
# add constraints on the anchor points
scalar = 1.2
if nadir[0] < 0:
m.const.add(expr=m.npv[0] >= scalar * nadir[0])
else:
m.const.add(expr=m.npv[0] >= 0 - scalar * nadir[0])
m.const.add(expr=m.total_LCA_midpoints[0, lca_type, lca_midpoint1] <= scalar * nadir[1])
m.const.add(expr=m.total_LCA_midpoints[0, lca_type, lca_midpoint2] <= scalar * nadir[2])
for i in range(3):
print(i, "utopia:", utopia[i], "nadir:", nadir[i])
# deactivate remaining objective function
m.Obj3.deactivate()
print(gwp, fe, npv)
return utopia, nadir, gwp, fe, npv
def aws(M, divs, midpoint, utopia, nadir, scenario, lca_type):
"""
the main control loop for conducting adaptive weight sums to identify the pareto front
:param M: the model being used
:param divs: the initial number of divisions
:param midpoint: the LCA ReCiPe lca_midpoint used
:param utopia: the utopia point for the pareto front
:param nadir: the nadir point for the pareto front
:param scenario: the scenario number
:param lca_type: the type of LCA being conducted (ALCA/CLCA)
:return: the npv values and the lca_midpoint values that comprise the pareto front
"""
# following Adaptive weighted-sum method for bi-objective optimization:Pareto front generation
print("\n\n AWS!!!!")
npv_new = []
gwp_new = []
# introduce the scaling factors
scalar = abs((nadir[0] - utopia[0]) / (nadir[1] - utopia[1]))
npv_vals, gwp_vals, models = ws(M, divs, lca_type, midpoint, scalar)
# append the anchor point to the pareto front
if len(npv_vals) > 0:
print("new pareto point appended - anchor", npv_vals[0], gwp_vals[0])
npv_new.append(npv_vals[0])
gwp_new.append(gwp_vals[0])
print_model(scenario, models[0], int(pyo.value(models[0].npv[0])), "TEA")
print_model(scenario, models[0], int(pyo.value(models[0].npv[0])), "LCA", midpoint=midpoint)
else:
# return if there are no values from weighted sums, there's no point in further refinement
return [], []
# calculate the euclidean distances between the points on the pareto front
dist = pareto_point_distance(gwp_vals, npv_vals)
# avoid those damn divide by zero errors
if len(dist) == 0:
return [], []
# refinement along sections of the pareto front
delta, n = ws_bound_refinement(dist, nadir, utopia)
# go through the points on the pareto front
for i in range(len(n)):
if n[i] > 1.6:
# determine the offset distances
try:
theta = math.atan(-((gwp_vals[i] - gwp_vals[i + 1]) / (npv_vals[i] - npv_vals[i + 1])))
except ZeroDivisionError:
# if we divide by zero, there is no need to subdivide this region any further
return [], []
delta1 = delta * math.cos(theta)
delta2 = delta * math.sin(theta)
print("\ndelta1", delta1)
print("delta2", delta2)
# solve the new submodel with additional constraints by calling the aws method
submodel = M.clone()
print("left point", npv_vals[i], gwp_vals[i])
print("right point", npv_vals[i + 1], gwp_vals[i + 1])
print("new lower npv bound", npv_vals[i] + delta1)
print("new upper gwp bound", gwp_vals[i + 1] + delta2)
submodel.const.add(expr=sum(submodel.npv[l] for l in submodel.Location) >= npv_vals[i] + delta1)
submodel.const.add(
expr=sum(submodel.total_LCA_midpoints[l, lca_type, midpoint] for l in submodel.Location) * scalar <=
gwp_vals[i + 1] + delta2)
# get results back from the submodel
x, y = aws(submodel, math.ceil(n[i]), midpoint, utopia, nadir, scenario, lca_type)
# append the returned values to our lists of points on the pareto front
for j in range(len(x)):
npv_new.append(x[j])
gwp_new.append(y[j])
# if the distance between this new point and the previous is less than half of delta, don't add the new point
pf_dist = (((npv_new[len(npv_new) - 1] - npv_vals[i + 1]) ** 2) +
((gwp_new[len(gwp_new) - 1] - gwp_vals[i + 1]) ** 2)) ** (1 / 2)
print("old point", npv_new[len(npv_new) - 1], gwp_new[len(gwp_new) - 1])
print("new point", npv_vals[i + 1], gwp_vals[i + 1])
print("distance along pareto front", pf_dist)
if pf_dist > delta / 2:
print("new pareto point appended!")
npv_new.append(npv_vals[i + 1])
gwp_new.append(gwp_vals[i + 1])
print_model(scenario, models[i + 1], int(pyo.value(models[i + 1].npv[0])), "TEA")
print_model(scenario, models[i + 1], int(pyo.value(models[i + 1].npv[0])), "LCA", midpoint=midpoint)
# because the pareto front is monotonic, we can sort the points without losing order (is it even necessary to sort?)
print("pareto front points", npv_new, gwp_new)
npv_new.sort()
gwp_new.sort()
print("pareto front points sorted", npv_new, gwp_new)
# return from the algorithm
print("\n\n\n\n\n return from an iteration!!!!!")
print(npv_new, gwp_new)
return npv_new, gwp_new
def pareto_point_distance(gwp_vals, npv_vals):
"""
calculates the distance between the pareto points
:param gwp_vals: the list of lca_midpoint values
:param npv_vals: the list of npv values
:return: the list continaing the distance between the points in the list
"""
dist = []
for i in range(len(npv_vals) - 1):
dist.append((((npv_vals[i] - npv_vals[i + 1]) ** 2) + ((gwp_vals[i] - gwp_vals[i + 1]) ** 2)) ** (1 / 2))
print("dist", dist)
return dist
def ws_bound_refinement(dist, nadir, utopia):
"""
Updates the lower and upper bounds on regions of the aws curve
:param dist: the list containing the distance between points
:param nadir: the nadir point for the pareto curve
:param utopia: the utopia point for the pareto curve
:return: the offset parameter delta and the number of refinements to conduct in each segment
"""
avg_length = sum(i for i in dist) / len(dist)
C = 1.5
delta = abs(nadir[0] - utopia[0]) / 8
if delta < 1000:
delta = 1000
print("delta", delta)
n = []
for i in range(len(dist)):
# to cut down on the number of iterations
if avg_length > delta:
temp = C * dist[i] / avg_length
else:
temp = C * dist[i] / delta
if temp > 6:
n.append(6)
else:
n.append(temp)
print("n", n)
return delta, n
def ws(M, divisions, lca_type, midpoint, scalar):
"""
Conducts weighted sums on a region of a pareto front
:param M: the model being optimized
:param divisions: the number of divisions to be made on the pareto front
:param lca_type: the LCA type (ALCA, CLCA)
:param midpoint: the ReCiPe lca_midpoint used on the pareto front
:param scalar: the scalar for the objective function so that the distance function works accurately
:return: a list of npv values, lca_midpoint values, and copies of the models for points alongside the pareto front
"""
# lists of objective values
models = []
npv_vals = []
gwp_vals = []
# do normal weighted sums
for i in range(divisions + 1):
alpha = 1 / divisions * i
print("\n\nalpha", alpha)
M.alpha = alpha
model = M
opt = pyo.SolverFactory('gurobi')
opt.options['TimeLimit'] = 600
try:
results = opt.solve(model, tee=True)
# check for infeasible solutions - don't add infeasible solutions to the pareto front
if not results.solver.termination_condition == TerminationCondition.optimal:
print("infeasible solution reached")
else:
npv_vals.append(pyo.value(model.npv[0]))
gwp_vals.append(pyo.value(model.total_LCA_midpoints[0, lca_type, midpoint]) * scalar)
models.append(model.clone())
except ValueError:
print("model did not find a solution within the time limit")
print("\n\n\nvalues from weighted sums iteration")
print("npv values", npv_vals)
print("gwp values", gwp_vals)
for i in range(len(models)):
print("model npv", i, pyo.value(models[i].npv[0]))
return npv_vals, gwp_vals, models
def ws3D(M, divisions, lca_type, midpoint1, midpoint2, s):
"""
Conducts weighted sums on a region of a pareto front
:param M: the model being optimized
:param divisions: the number of divisions to be made on the pareto front
:param lca_type: the LCA type (ALCA, CLCA)
:param midpoint1: the ReCiPe lca_midpoint used on the pareto front
:param midpoint2: the ReCiPe lca_midpoint used on the pareto front
:param s: the scenario number for the problem
:return: a list of npv values, lca_midpoint values, and copies of the models for points alongside the pareto front
"""
# lists of objective values
x = []
y = []
z = []
# do normal weighted sums
for i in range(divisions + 1):
# update objective function parameters
alpha = 1 / divisions * i
M.alpha = alpha
for j in range(divisions + 1):
# update objective function parameters
alpha2 = 1 / divisions * j
print("\n\nalpha", alpha)
print("alpha2", alpha2)
M.alpha2 = alpha2
# solve the model
model = M
opt = pyo.SolverFactory('gurobi')
opt.options['TimeLimit'] = 600
opt.options["mipgap"] = .0008
try:
results = opt.solve(model, tee=True)
# check for infeasible solutions - don't add infeasible solutions to the pareto front
if not results.solver.termination_condition == TerminationCondition.optimal:
print("infeasible solution reached")
else:
# save data for plotting
x.append(pyo.value(model.npv[0]))
y.append(pyo.value(model.total_LCA_midpoints[0, lca_type, midpoint1]))
z.append(pyo.value(model.total_LCA_midpoints[0, lca_type, midpoint2]))
# write out excel data to reduce memory usage
print_model(s, model, int(pyo.value(model.npv[0])), "TEA")
print_model(s, model, int(pyo.value(model.npv[0])), "LCA", midpoint=midpoint1)
except ValueError:
print("model did not find a solution within the time limit")
# log outputs
print("\n\n\nvalues from weighted sums iteration")
print("npv values", x)
print("gwp values", y)
print("fe values", z)
return x, y, z
def pareto_front2D(M, midpoint, scenario, A, lca_type):
"""
Calculates the pareto front
:param M: the model used in the pareto front
:param midpoint: the lca_midpoint used for the objective function
:param scenario: the scenario number
:param A: parameters from the CLCO_Data
:param lca_type: the LCA type (ALCA/CLCA)
:return:
"""
M.Obj = pyo.Objective(expr=sum(M.npv[l] for l in M.Location), sense=pyo.maximize)
M.Obj2 = pyo.Objective(expr=sum(M.total_LCA_midpoints[l, lca_type, midpoint] for l in M.Location),
sense=pyo.minimize)
try:
utopia, nadir = utopian(M, midpoint, lca_type)
except ValueError as err:
print(err.args)
print("no further attempt to find the pareto front")
return
print("\n\n finished computing utopia and nadir points")
# use the new objective function with the new weights
M.combined = pyo.Objective(
expr=M.alpha * sum(M.npv[l] for l in M.Location) -
(1 - M.alpha) * abs((nadir[0] - utopia[0]) / (nadir[1] - utopia[1])) * sum(
M.total_LCA_midpoints[l, lca_type, midpoint] for l in M.Location), sense=pyo.maximize)
print("returned to control method")
# gather the points on the pareto front
x, y = aws(M, 4, midpoint, utopia, nadir, scenario, lca_type)
# rescale the y points back to their original values
scalar = abs((nadir[0] - utopia[0]) / (nadir[1] - utopia[1]))
y_rescaled = [sum(i / (scalar * A.FEEDSTOCK_SUPPLY[l] * A.TIME_PERIODS) for l in M.Location) for i in y]
x_rescaled = [sum(i / (A.FEEDSTOCK_SUPPLY[l] * A.TIME_PERIODS) for l in M.Location) for i in x]
print("rescaled points")
# plot the points
plt.clf()
plt.plot(x_rescaled, y_rescaled, 'ob-')
plt.title("Pareto Front (Adaptive Weighted Sums)")
plt.xlabel("NPV ($USD) per ton manure")
plt.ylabel(str(midpoint) + " impact per ton manure")
plt.savefig(save_plot(scenario, midpoint=midpoint), dpi=300)
print("saved fig")
# plt.show()
return 1
def pareto_front3D(M, midpoint1, midpoint2, scenario, A, lca_type):
"""
Calculates the 3D pareto front
:param M: the model used in the pareto front
:param midpoint1: the first lca_midpoint used for the objective function
:param midpoint2: the second lca_midpoint used for the objective function
:param scenario: the scenario number
:param A: parameters from the CLCO_Data
:param lca_type: the LCA type (ALCA/CLCA)
:return:
"""
print("3D pareto front!!!!")
M.Obj = pyo.Objective(expr=0 - sum(M.npv[l] for l in M.Location), sense=pyo.minimize)
M.Obj2 = pyo.Objective(expr=sum(M.total_LCA_midpoints[l, lca_type, midpoint1] for l in M.Location),
sense=pyo.minimize)
M.Obj3 = pyo.Objective(expr=sum(M.total_LCA_midpoints[l, lca_type, midpoint2] for l in M.Location),
sense=pyo.minimize)
try:
utopia, nadir, g, f, n = utopian3D(M, midpoint1, midpoint2, lca_type)
except ValueError as err:
print(err.args)
print("no further attempt to find the pareto front")
return
print("\n\n finished computing utopia and nadir points")
# use the new objective function with the new weights
ranges = [abs(nadir[0] - utopia[0]), abs(nadir[1] - utopia[1]), abs(nadir[2] - utopia[2])]
print(ranges)
M.combined = pyo.Objective(
expr=0 - M.alpha * M.alpha2 * (sum(M.npv[l] for l in M.Location) * max(ranges) / ranges[0]) +
(1 - M.alpha) * M.alpha2 * (
sum(M.total_LCA_midpoints[l, lca_type, midpoint1] for l in M.Location) * max(ranges) / ranges[1]) +
(1 - M.alpha2) * (
sum(M.total_LCA_midpoints[l, lca_type, midpoint2] for l in M.Location) * max(ranges) / ranges[2])
, sense=pyo.minimize)
print("returned to control method")
# gather the points on the pareto front
x, y, z = ws3D(M, 6, lca_type, midpoint1, midpoint2, scenario)
for j in range(len(g)):
x.append(n[j])
y.append(g[j])
z.append(f[j])
# x, y, z = GPBAB(M, 7, midpoint1, midpoint2, lca_type, ranges, utopia, nadir)
# rescale the y points back to their original values
npv_rescaled = [i / (A.FEEDSTOCK_SUPPLY[0] * A.TIME_PERIODS) for i in x]
gwp_rescaled = [i / (A.FEEDSTOCK_SUPPLY[0] * A.TIME_PERIODS) for i in y]
fe_rescaled = [i / (A.FEEDSTOCK_SUPPLY[0] * A.TIME_PERIODS) for i in z]
print("rescaled points")
# plot the points
plt.clf()
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
surf = ax.plot_trisurf(gwp_rescaled, fe_rescaled, npv_rescaled, linewidth=0.1, cmap=plt.cm.CMRmap)
fig.colorbar(surf, shrink=0.5, aspect=5)
ax.set_zlabel("NPV ($USD) per ton manure")
ax.set_xlabel("GWP (kg CO2-eq) per ton manure")
ax.set_ylabel("FE (kg P-eq) per ton manure")
# plt.show()
plt.clf()
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.stem(gwp_rescaled, fe_rescaled, npv_rescaled, markerfmt=" ", basefmt=" ", linefmt="grey",
bottom=min(npv_rescaled) - 3)
sc = ax.scatter(gwp_rescaled, fe_rescaled, npv_rescaled, s=80, cmap='plasma', c=npv_rescaled, edgecolor="black",
alpha=1)
fig.colorbar(sc, shrink=0.5, aspect=5)
ax.set_zlabel("NPV ($USD) per ton manure")
ax.set_xlabel("GWP (kg CO2-eq) per ton manure")
ax.set_ylabel("FE (kg P-eq) per ton manure")
# plt.show()
return 1
def GPBAB(M, delta, midpoint1, midpoint2, lca_type, r, u, n):
"""
alternative pareto front generation algorithm - untested
:param M: the model being used
:param delta: the number of points to be spread across the solution space
:param midpoint1: the first lca midpoint to use
:param midpoint2: the second lca midpoint to use
:param lca_type: the lca type being used for the optimization algorithm
:param r: the range of parameter values
:param u: the utopia points for each objective
:param n: the nadir points for each objective
:return: none - program is still untested
"""
# initialize loop control variables
e2 = n[2]
offset1 = r[1] / delta
offset2 = r[2] / delta
zwv2 = u[2]
M.alpha = 0.999
M.alpha2 = 0.999
x = []
y = []
z = []
j = 0
while e2 > u[2]:
i = 0
j = j + 1
e1 = n[1]
while e1 > u[1]:
i = i + 1
# add constraints
M.const.add(expr=sum(M.total_LCA_midpoints[l, lca_type, midpoint1] for l in M.Location) <= e1)
M.const.add(expr=sum(M.total_LCA_midpoints[l, lca_type, midpoint2] for l in M.Location) <= e2)
# update the epsilon bounds
print("\n\ne1", e1)
print("e2", e2)
print("zwv2", zwv2)
model = M
opt = pyo.SolverFactory('gurobi')
opt.options['TimeLimit'] = 180
try:
results = opt.solve(model, tee=True)
# check for infeasible solutions - don't add infeasible solutions to the pareto front
if not results.solver.termination_condition == TerminationCondition.optimal:
print("infeasible solution reached")
else:
# get results from model
x.append(pyo.value(model.npv[0]))
y.append(pyo.value(model.total_LCA_midpoints[0, lca_type, midpoint1]))
z.append(pyo.value(model.total_LCA_midpoints[0, lca_type, midpoint2]))
print_model(scenario, model, i + int(pyo.value(model.npv[0])), "TEA")
print_model(scenario, model, i + int(pyo.value(model.npv[0])), "LCA", midpoint=midpoint1)
# update parameters
e1 = pyo.value(model.total_LCA_midpoints[0, lca_type, midpoint1]) - offset1
zwv2 = max(zwv2, pyo.value(model.total_LCA_midpoints[0, lca_type, midpoint2]))
except ValueError:
print("model did not find a solution within the time limit")
# therefore, any future refinements are futile
e1 = u[1] - 1
# update the outer loop
e2 = zwv2 - offset2
def initialize_model(scenario, j, midpoint, lca_type):
"""
Initializes the CLCO model
:param scenario: the scenario number
:param j: the county number
:param midpoint: the lca_midpoint type for MOO
:param lca_type: the LCA type for MOO
:return:
"""
A = CLCO_Data(scenario)
M = pyo.ConcreteModel(scenario)
add_sets(A, M, j)
add_variables(M)
add_constraints(A, M, scenario)
# solving the model
if 1000 < scenario < 3000:
return pareto_front2D(M, midpoint, scenario, A, lca_type)
elif scenario in [4501, 4502, 4503, 4511, 4512, 4513, 4521, 4522, 4523]:
return pareto_front3D(M, 'climate change', 'eutrophication: freshwater', scenario, A, lca_type)
elif 2999 < scenario < 9999:
if int((scenario / 100) % 10) == 0:
M.Obj = pyo.Objective(expr=sum(M.npv[l] for l in M.Location), sense=pyo.maximize)
if int((scenario / 100) % 10) == 1:
M.Obj = pyo.Objective(
expr=sum(M.npv[l] - M.total_LCA_midpoints[l, "CLCA", "climate change"] / 5 for l in M.Location),
sense=pyo.maximize)
if int((scenario / 100) % 10) == 2:
M.Obj = pyo.Objective(
expr=sum(M.npv[l] - M.total_LCA_midpoints[l, "CLCA", "climate change"] for l in M.Location),
sense=pyo.maximize)
elif scenario == 50:
M.Obj = pyo.Objective(expr=sum(M.total_LCA_midpoints[l, lca_type, "climate change"] for l in M.Location),
sense=pyo.minimize)
elif scenario in [425]:
M.Obj = pyo.Objective(
expr=sum(M.npv[l] - M.total_LCA_midpoints[l, "CLCA", "climate change"] / 5 for l in M.Location),
sense=pyo.maximize)
else:
M.Obj = pyo.Objective(expr=sum(M.npv[l] for l in M.Location), sense=pyo.maximize)
model = M
opt = pyo.SolverFactory('gurobi')
if scenario in [5]:
opt.options["mipgap"] = .001
print(opt.solve(model, tee=True))
model.npv.pprint()
# print data to csv
print_model(scenario, model, j, "TEA")
print_model(scenario, model, j, "LCA")
def add_constraints(A, M, scenario):
"""
Adds constraints to the model
:param A: parameters for the model
:param M: the model
:param scenario: the scenario number
:return:
"""
# constraint list definition
M.const = pyo.ConstraintList()
for l in M.Location:
# only one stage can be selected for each location for AD
M.const.add(expr=sum(M.decision_ad_stage[l, stage] for stage in M.ADStages) == 1)
# time dependent constraints
for t in M.Time:
facility_constraints(A, M, l, scenario, t)
revenue_constraints(A, M, l, scenario, t)
opex_constraints(A, M, l, t)
# time independent constraints
capex_constraints(A, M, l)
lca_constraints(A, M, l)
# NPV discounting for each location
M.const.add(expr=M.npv[l] == - sum(M.process_capex[l, tech] + M.storage_capex[l, tech] for tech in M.Technology)
+ sum(M.opex_revenues[l, t, tech, revenue_type]
for t in M.Time for tech in M.Technology for revenue_type in M.OPEXSubRevenues) -
sum(M.opex_costs[l, t, tech, cost_type]
for cost_type in M.OPEXSubCosts for t in M.Time for tech in M.Technology))
def lca_constraints(A, M, l):
"""
Adds constraints on ALCA/CLCA to model
:param A: data
:param M: model
:param l: location
:return: N/A
"""
# LCA calculations
for cat in M.LCAMidpointCat:
for lca_type in M.LCATypes:
# mid point for a point in time
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'natural gas', cat] ==
A.IMPACT[lca_type, 'natural gas', cat] * sum(M.purchased_fuel[l, t, tech]
for tech in M.Technology for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'grid electricity', cat] ==
A.IMPACT[lca_type, 'grid electricity', cat] * sum(M.purchased_power[l, t, tech]
for tech in M.Technology for t in
M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'diesel', cat] == A.IMPACT[lca_type, 'diesel', cat]
* (sum(M.inputs[l, t, tech, 'bio-oil diesel'] for tech in M.Technology for t in M.Time) +
sum(M.inputs[l, t, tech, 'transportation'] * A.INTRA_COUNTY_TRANSPORT_DISTANCE[l]
* A.DIESEL_USE for t in M.Time for tech in M.Technology)))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'water', cat] == A.IMPACT[lca_type, 'water', cat] *
sum(M.inputs[l, t, tech, 'water'] for tech in M.Technology for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'biochar-chp', cat] ==
sum(A.IMPACT[lca_type, 'biochar-chp', temp, cat] * M.biochar_from_pyrolysis[
l, t, feed, temp, 'CHP'] for temp in M.PyrolysisTemperatures
for feed in M.PyrolysisFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'biochar-land', cat] ==
sum(A.IMPACT[lca_type, 'biochar-land', temp, cat] * M.biochar_from_pyrolysis[
l, t, feed, temp, 'land'] for temp in M.PyrolysisTemperatures
for feed in M.PyrolysisFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'biochar-disposal', cat] ==
A.IMPACT[lca_type, 'biochar-disposal', cat] *
sum(M.biochar_from_pyrolysis[l, t, feed, temp, 'disposal']
for temp in M.PyrolysisTemperatures for feed in M.PyrolysisFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'pyro-bio-oil-chp', cat] ==
sum(A.IMPACT[lca_type, 'pyro-bio-oil-chp', temp, cat] * M.biooil_from_pyrolysis[
l, t, feed, temp, 'CHP'] for temp in M.PyrolysisTemperatures for feed
in M.PyrolysisFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'syngas-chp', cat] ==
sum(A.IMPACT[lca_type, 'syngas-chp', temp, cat] * M.syngas_from_pyrolysis[
l, t, feed, temp, 'CHP'] for temp in M.PyrolysisTemperatures
for feed in M.PyrolysisFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'syngas-disposal', cat] ==
sum(A.IMPACT[lca_type, 'syngas-disposal', temp, cat] * M.syngas_from_pyrolysis[
l, t, feed, temp, 'disposal'] for temp in M.PyrolysisTemperatures
for feed in M.PyrolysisFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'pyro-ap-disposal', cat] ==
A.IMPACT[lca_type, 'pyro-ap-disposal', cat] *
sum(M.ap_from_pyrolysis[l, t, feed, temp, 'disposal']
for temp in M.PyrolysisTemperatures for feed in M.PyrolysisFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htl-hydrochar-land', cat] ==
A.IMPACT[lca_type, 'htl-hydrochar-land', cat]
* sum(M.hydrochar_from_htl[l, t, feed, temp, 'land']
for temp in M.HTLTemperatures for feed in M.HTLFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htl-hydrochar-chp', cat] ==
A.IMPACT[lca_type, 'htl-hydrochar-chp', cat]
* sum(M.hydrochar_from_htl[l, t, feed, temp, 'CHP']
for temp in M.HTLTemperatures for feed in M.HTLFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htl-hydrochar-disposal', cat] ==
A.IMPACT[lca_type, 'htl-hydrochar-disposal', cat] *
sum(M.hydrochar_from_htl[l, t, feed, temp, 'disposal'] for temp in M.HTLTemperatures
for feed in M.HTLFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htl-bio-oil-chp', cat] ==
A.IMPACT[lca_type, 'htl-bio-oil-chp', cat] * sum(M.biooil_from_htl[l, t, feed, temp, 'CHP']
for temp in M.HTLTemperatures for feed in
M.HTLFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htl-gp-disposal', cat] ==
A.IMPACT[lca_type, 'htl-gp-disposal', cat] *
sum(M.gp_from_htl[l, t, feed, temp, 'disposal']
for temp in M.HTLTemperatures for feed in M.HTLFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htl-ap-disposal', cat] ==
A.IMPACT[lca_type, 'htl-ap-disposal', cat] *
sum(M.ap_from_htl[l, t, feed, temp, 'disposal']
for temp in M.HTLTemperatures for feed in M.HTLFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htc-hydrochar-land', cat] ==
A.IMPACT[lca_type, 'htc-hydrochar-land', cat] * sum(
M.hydrochar_from_htc[l, t, feed, temp, 'land']
for temp in M.HTCTemperatures for feed in M.HTCFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htc-hydrochar-chp', cat] ==
sum(A.IMPACT[lca_type, 'htc-hydrochar-chp', temp, cat] * M.hydrochar_from_htc[
l, t, feed, temp, 'CHP']
for temp in M.HTCTemperatures for feed in M.HTCFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htc-hydrochar-disposal', cat] ==
A.IMPACT[lca_type, 'htc-hydrochar-disposal', cat] *
sum(M.hydrochar_from_htc[l, t, feed, temp, 'disposal']
for temp in M.HTCTemperatures for feed in M.HTCFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htc-gp-disposal', cat] ==
sum(A.IMPACT[lca_type, 'htc-gp-disposal', temp, cat] *
M.gp_from_htc[l, t, feed, temp, 'disposal']
for temp in M.HTCTemperatures for feed in M.HTCFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'htc-ap-disposal', cat] ==
A.IMPACT[lca_type, 'htc-ap-disposal', cat] *
sum(M.ap_from_htc[l, t, feed, temp, 'disposal']
for temp in M.HTCTemperatures for feed in M.HTCFeedstocks for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'digestate-land', cat] ==
A.IMPACT[lca_type, 'digestate-land', cat] * sum(M.digestate_from_ad[l, t, stage, 'land']
for stage in M.ADStages for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'digestate-disposal', cat] ==
A.IMPACT[lca_type, 'digestate-disposal', cat] *
sum(M.digestate_from_ad[l, t, stage, 'disposal'] for stage in M.ADStages for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'biogas-disposal', cat] ==
A.IMPACT[lca_type, 'biogas-disposal', cat] * sum(M.biogas_from_ad[l, t, stage, 'disposal']
for stage in M.ADStages for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'biogas-chp', cat] ==
A.IMPACT[lca_type, 'biogas-chp', cat] * sum(M.biogas_from_ad[l, t, stage, 'CHP']
for stage in M.ADStages for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'manure-land', cat] ==
A.IMPACT[lca_type, 'manure-land', cat] *
sum(M.feedstock_from_storage[l, t] for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'facility construction', cat] ==
A.IMPACT[lca_type, 'facility construction', "ad", cat] * M.process_capacity[l, 'AD'] +
A.IMPACT[lca_type, 'facility construction', "chem", cat] * (M.process_capex[l, 'Pyrolysis']
+ M.process_capex[l, 'HTL'] +
M.process_capex[l, 'HTC']))
# CHP facility construction LCA impacts are spread out across the various CHP categories
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'storage-facility-solids', cat] ==
A.IMPACT[lca_type, 'storage-facility-solids', cat] * (M.feedstock_storage_capacity[l] +
M.pyrolysis_storage_capacity[
l, 'Biochar'] +
M.htl_storage_capacity[
l, 'Hydrochar'] +
M.htc_storage_capacity[
l, 'Hydrochar'] +
M.ad_storage_capacity[
l, 'digestate']))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'storage-facility-liquids', cat] ==
A.IMPACT[lca_type, 'storage-facility-liquids', cat] *
(M.pyrolysis_storage_capacity[l, 'AP'] +
M.pyrolysis_storage_capacity[l, 'Syngas'] +
M.htl_storage_capacity[l, 'AP'] +
M.htc_storage_capacity[l, 'AP'] +
M.ad_storage_capacity[l, 'biogas']))
M.const.add(expr=M.LCA_midpoints[l, lca_type, "biochar market", cat] == 0)
'''A.IMPACT[lca_type,
'biochar market', cat] *
sum(M.biochar_from_pyrolysis[l, t, feedstock, temp, 'market']
for l in M.Location for t in M.Time for feedstock in M.PyrolysisFeedstocks
for temp in M.PyrolysisTemperatures))''' # for when biochar can be sold on the market - still needs to be implemented
M.const.add(expr=M.LCA_midpoints[l, lca_type, "bio-oil market", cat] ==
A.IMPACT[lca_type, 'bio-oil market', cat] *
(sum(M.biooil_from_pyrolysis[l, t, feedstock, temp, 'market']
for l in M.Location for t in M.Time for feedstock in M.PyrolysisFeedstocks
for temp in M.PyrolysisTemperatures)
+ sum(M.biooil_from_htl[l, t, feedstock, temp, 'market']
for l in M.Location for t in M.Time for feedstock in M.HTLFeedstocks
for temp in M.HTLTemperatures)))
M.const.add(expr=M.LCA_midpoints[l, lca_type, "hydrochar market", cat] ==
A.IMPACT[lca_type, 'hydrochar market', cat] *
(sum(M.hydrochar_from_htc[l, t, feedstock, temp, 'market']
for l in M.Location for t in M.Time for feedstock in M.HTCFeedstocks
for temp in M.HTCTemperatures)
+ sum(M.hydrochar_from_htl[l, t, feedstock, temp, 'market']
for l in M.Location for t in M.Time for feedstock in M.HTLFeedstocks
for temp in M.HTLTemperatures)))
M.const.add(expr=M.total_LCA_midpoints[l, lca_type, cat] ==
sum(M.LCA_midpoints[l, lca_type, origin, cat] for origin in M.ALCAInputs))
# LCA specific constraints
if lca_type == "ALCA":
M.const.add(expr=M.LCA_midpoints[l, lca_type, "avoided electricity", cat] == 0)
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'N fertilizer', cat] == 0)
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'P fertilizer', cat] == 0)
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'K fertilizer', cat] == 0)
elif lca_type == "CLCA":
M.const.add(expr=M.LCA_midpoints[l, lca_type, "avoided electricity", cat] ==
-A.IMPACT[lca_type, 'grid electricity', cat] *
sum(M.chp_market[l, t, 'electricity'] for t in M.Time))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'N fertilizer', cat] ==
-A.IMPACT[lca_type, 'N fertilizer', cat] *
sum(M.avoided_fertilizers[l, t, tech, 'N']
for l in M.Location for t in M.Time for tech in M.Technology))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'P fertilizer', cat] ==
-A.IMPACT[lca_type, 'P fertilizer', cat] *
sum(M.avoided_fertilizers[l, t, tech, 'P']
for l in M.Location for t in M.Time for tech in M.Technology))
M.const.add(expr=M.LCA_midpoints[l, lca_type, 'K fertilizer', cat] ==
-A.IMPACT[lca_type, 'K fertilizer', cat] *
sum(M.avoided_fertilizers[l, t, tech, 'K']
for l in M.Location for t in M.Time for tech in M.Technology))
def capex_constraints(A, M, l):
"""
Constraints for calculating the CAPEX for the model
:param A: parameter information
:param M: the model
:param l: the location
:return: N/A
"""
# CAPEX for storage and process capacity
# implementing piecewise linear approximation
q = symbols("q")
lpa_xvals = []
[lpa_xvals.append(x) for x in A.ORIGINAL_FEEDSTOCK_SUPPLY if x not in lpa_xvals]
lpa_xvals.append(0.00000001)
for i in range(16):
lpa_xvals.append(i * 1000 + 4500)
for i in range(19):
lpa_xvals.append(i * 10000 + 20000)
lpa_xvals.sort()
pyro_process = A.CAPEX['Pyrolysis', 'process', 'coefficient'] * (q * 1000 / A.HOURS_PER_PERIOD) ** A.CAPEX[
'Pyrolysis', 'process', 'exponent']
pyro_storage = A.CAPEX['Pyrolysis', 'storage', 'coefficient'] * (q / A.HOURS_PER_PERIOD) ** A.CAPEX[
'Pyrolysis', 'storage', 'exponent']
htl_process = A.CAPEX['HTL', 'process', 'coefficient'] * (q / A.DRY_BIOMASS_REF) ** A.CAPEX[
'HTL', 'process', 'exponent']
htl_storage = A.CAPEX['HTL', 'storage', 'coefficient'] * (q / A.HOURS_PER_PERIOD) ** A.CAPEX[
'HTL', 'storage', 'exponent']
htc_process = A.CAPEX['HTC', 'process', 'coefficient'] * (q / A.DRY_BIOMASS_REF) ** A.CAPEX[
'HTC', 'process', 'exponent']
htc_storage = A.CAPEX['HTC', 'storage', 'coefficient'] * (q / A.HOURS_PER_PERIOD) ** A.CAPEX[
'HTC', 'storage', 'exponent']
ad_process = A.CAPEX['AD', 'process', 'coefficient'] * q ** A.CAPEX['AD', 'process', 'exponent']
ad_storage = A.CAPEX['AD', 'storage', 'coefficient'] * q ** A.CAPEX['AD', 'storage', 'exponent']
chp_process = A.CAPEX['CHP', 'process', 'coefficient'] * q ** A.CAPEX['CHP', 'process', 'exponent']
solid_storage = A.CAPEX['Solid', 'storage', 'coefficient'] * q ** A.CAPEX['Solid', 'storage', 'exponent']
pyro_proc_lpa = [pyro_process.evalf(subs={q: x}) for x in lpa_xvals]
pyro_stor_lpa = [pyro_storage.evalf(subs={q: x}) for x in lpa_xvals]
solids_products = [solid_storage.evalf(subs={q: x}) for x in lpa_xvals]
htl_proc_lpa = [htl_process.evalf(subs={q: x}) for x in lpa_xvals]
htl_stor_lpa = [htl_storage.evalf(subs={q: x}) for x in lpa_xvals]
htc_proc_lpa = [htc_process.evalf(subs={q: x}) for x in lpa_xvals]
htc_stor_lpa = [htc_storage.evalf(subs={q: x}) for x in lpa_xvals]
ad_proc_lpa = [ad_process.evalf(subs={q: x}) for x in lpa_xvals]
ad_stor_lpa = [ad_storage.evalf(subs={q: x}) for x in lpa_xvals]
chp_proc_lpa = [chp_process.evalf(subs={q: x}) for x in lpa_xvals]
# Piecewise linear approximation to capital cost terms
for tech in M.Technology:
if tech == 'Pyrolysis':
M.con5 = pyo.Piecewise(M.process_capex[l, 'Pyrolysis'], M.process_capacity[l, 'Pyrolysis'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=pyro_proc_lpa,
pw_repn='SOS2')
M.con6 = pyo.Piecewise(M.pyro_storage_cost[l, 'Biochar'], M.pyrolysis_storage_capacity[l, 'Biochar'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=solids_products,
pw_repn='SOS2')
M.con7 = pyo.Piecewise(M.pyro_storage_cost[l, 'Biooil'], M.pyrolysis_storage_capacity[l, 'Biooil'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=pyro_stor_lpa,
pw_repn='SOS2')
M.con8 = pyo.Piecewise(M.pyro_storage_cost[l, 'AP'], M.pyrolysis_storage_capacity[l, 'AP'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=pyro_stor_lpa,
pw_repn='SOS2')
M.con9 = pyo.Piecewise(M.pyro_storage_cost[l, 'Syngas'], M.pyrolysis_storage_capacity[l, 'Syngas'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=pyro_stor_lpa,
pw_repn='SOS2')
M.const.add(expr=M.storage_capex[l, 'Pyrolysis'] == sum(
M.pyro_storage_cost[l, prod] for prod in M.PyrolysisProducts))
elif tech == 'HTL':
M.con10 = pyo.Piecewise(M.process_capex[l, 'HTL'], M.process_capacity[l, 'HTL'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=htl_proc_lpa,
pw_repn='SOS2')
M.con11 = pyo.Piecewise(M.htl_storage_cost[l, 'Hydrochar'], M.htl_storage_capacity[l, 'Hydrochar'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=solids_products,
pw_repn='SOS2')
M.con12 = pyo.Piecewise(M.htl_storage_cost[l, 'Biooil'], M.htl_storage_capacity[l, 'Biooil'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=htl_stor_lpa,
pw_repn='SOS2')
M.con13 = pyo.Piecewise(M.htl_storage_cost[l, 'AP'], M.htl_storage_capacity[l, 'AP'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=htl_stor_lpa,
pw_repn='SOS2')
M.con14 = pyo.Piecewise(M.htl_storage_cost[l, 'GP'], M.htl_storage_capacity[l, 'GP'],
pw_pts=lpa_xvals,
pw_constr_type='EQ',
f_rule=htl_stor_lpa,
pw_repn='SOS2')
M.const.add(expr=M.storage_capex[l, 'HTL'] == sum(M.htl_storage_cost[l, prod] for prod in M.HTLProducts))