-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathrhp_plan.py
118 lines (101 loc) · 6.86 KB
/
rhp_plan.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
import numpy as np
from copy import deepcopy
from scipy import optimize
#See the rhp_intro.pdf for an explanation of the comments
#Plan function; full designates including exports
def plan(time_steps, planning_horizon, primary_resource_list, supply_use_list, use_imported_list,
depreciation_list, full_domestic_target_output_list, imported_target_output_list,
export_constraint_boolean, export_prices_list, import_prices_list, upper_bound_on_activity,
max_iterations, tolerance):
result_list, lagrange_list, slack_list = [], [], []
aggregate_constraint_matrix_final_list = []
aggregate_constraint_vector_final_list = []
aggregate_primary_resource_final_list = []
modifier_value = deepcopy(np.zeros_like(full_domestic_target_output_list[0]))
if export_constraint_boolean:
#Import prices list is a list of vectors, but it is converted to matrices where the
# diagonal takes on the vectors values and the other elements are 0
import_prices_list = deepcopy([np.diagflat(import_prices_list[i]) for i in range(len(import_prices_list))])
for N in range(time_steps):
#Constructing the aggregate constraint matrix (each DJ)
zero_matrix = np.zeros_like(supply_use_list[0])
horizontal_block_list = []
supply_use_part = deepcopy(supply_use_list[N])
for i in range(planning_horizon+1):
if (i < planning_horizon):
zero_part = deepcopy(np.block([zero_matrix]*(planning_horizon-i)))
horizontal_block_list.append(np.block([supply_use_part, zero_part]))
#depreciation
supply_use_part = deepcopy(np.matmul(depreciation_list[N+i],supply_use_part))
supply_use_part = deepcopy(np.block([supply_use_part,supply_use_list[N+i]]))
horizontal_block_list.append(supply_use_part)
aggregate_constraint_matrix = np.vstack(horizontal_block_list)
#Constructing the aggregate constraint vector (each Dr)
vertical_block_list = []
modifier = deepcopy(-modifier_value)
for i in range(planning_horizon+1):
vertical_block_list.append(np.sum([full_domestic_target_output_list[N+i], modifier], axis=0))
#depreciation and carry
modifier = deepcopy(np.matmul(depreciation_list[N+i], vertical_block_list[-1]))
aggregate_constraint_vector = np.vstack(vertical_block_list)
#Constructing the aggregate primary resource vector (each c)
aggregate_primary_resource_list = []
for i in range(planning_horizon+1):
aggregate_primary_resource_list.append(primary_resource_list[N+i])
aggregate_primary_resource_vector = np.vstack(aggregate_primary_resource_list)
#Export constraint
if export_constraint_boolean:
#Constructing the use imports constraint matrix (each p_imp^T T-bar)
horizontal_block_list = []
for i in range(planning_horizon+1):
price_augmented_use_imports = deepcopy(np.matmul(import_prices_list[N+i],
np.transpose(use_imported_list[N+i])))
horizontal_block_list.append(np.transpose(-price_augmented_use_imports))
use_imports_constraint_matrix = np.hstack(horizontal_block_list)
#Constructing the export constraint matrix (each p_exp^T)
horizontal_block_list = []
for i in range(planning_horizon+1):
export_price_diagonal = deepcopy(np.transpose(export_prices_list[N+i]))
horizontal_block_list.append(np.diagflat(export_price_diagonal))
aggregate_export_price_matrix = np.hstack(horizontal_block_list)
zero_like = np.zeros_like(aggregate_export_price_matrix)
aggregate_zero_matrix = np.tile(zero_like, (planning_horizon+1,1))
export_constraint_matrix = np.vstack([aggregate_zero_matrix, aggregate_export_price_matrix])
#Constructing the directly imported target output constraint vector (p_imp^T r_imp)
vertical_block_list = []
for i in range(planning_horizon+1):
vertical_block_list.append(np.matmul(import_prices_list[N+i], imported_target_output_list[N+i]))
imported_target_output_constraint_vector = np.sum(vertical_block_list, axis=0)
#Constructing the export augmentet constraint objects
aggregate_constraint_matrix = deepcopy(np.vstack([aggregate_constraint_matrix, use_imports_constraint_matrix]))
aggregate_constraint_matrix = deepcopy(np.hstack([aggregate_constraint_matrix, export_constraint_matrix]))
aggregate_constraint_vector = deepcopy(np.vstack([aggregate_constraint_vector, imported_target_output_constraint_vector]))
ones_length = aggregate_primary_resource_vector.shape[0]
one_vector = deepcopy(np.matrix(np.ones((ones_length, 1))))
aggregate_primary_resource_vector = deepcopy(np.vstack([aggregate_primary_resource_vector, one_vector]))
# Plan
result = optimize.linprog(c=aggregate_primary_resource_vector,
A_ub=-aggregate_constraint_matrix,
b_ub=-aggregate_constraint_vector, options = {'maxiter': max_iterations, 'tol': tolerance},
bounds=(0, upper_bound_on_activity), method='highs')
print(result.success)
print(result.status)
lagrange_ineq = -optimize.linprog(c=aggregate_primary_resource_vector,
A_ub=-aggregate_constraint_matrix,
b_ub=-aggregate_constraint_vector,
bounds=(0, upper_bound_on_activity), options = {'maxiter': 1000, 'tol': 1e-8},
method='highs')['ineqlin']['marginals']
result_list.append(result.x)
lagrange_list.append(lagrange_ineq)
slack_list.append(result.slack)
aggregate_constraint_matrix_final_list.append(aggregate_constraint_matrix)
aggregate_constraint_vector_final_list.append(aggregate_constraint_vector)
aggregate_primary_resource_final_list.append(aggregate_primary_resource_vector)
#Production carry into the next time step
if export_constraint_boolean:
recent_slack = deepcopy(np.array_split(slack_list[N], (planning_horizon+2)))
modifier_value = deepcopy(np.matmul(depreciation_list[N], recent_slack[0]).reshape([-1,1]))
else:
recent_slack = deepcopy(np.array_split(slack_list[N], planning_horizon+1))
modifier_value = deepcopy(np.matmul(depreciation_list[N], recent_slack[0]).reshape([-1,1]))
return([result_list, lagrange_list, slack_list, aggregate_primary_resource_final_list, aggregate_constraint_matrix_final_list, aggregate_constraint_vector_final_list])