Это мой код. Я пытаюсь оптимизировать количество воды для распределения от узлов предложения до спроса по транспортной модели. Есть два ограничения: спрос должен быть как минимум обеспечен, и вы не можете взять больше воды, чем тот, который он способен распределить; но модель не соблюдает это ограничение.
import numpy as np
from scipy.optimize import minimize, LinearConstraint, Bounds
# Data
rho = 1000
g = 9.8
C = 140.0 #roughness coefficient
J = rho*g/(float(1e6)) #constant of function
capacity = [0.3875, 0.607, 0.374, 2.404, 3.614, 2.738, 0.037, 0.17, 0.17, 2.63] #Supply data
demand = [0.768, 0.315, 0.38, 0.463, 0.094, 0.46917, 0.095, 0.129, 0.24057, 0.55, 0.03981, 0.518, 0.5, 0.5, 0.13, 0.021, 0.03721, 0.01346] #Demand data
m= len(capacity)
n = len(demand)
x0 = np.array(m*n*[1]) #Initial point for algorithm
# In[59]:
#READ L and h FROM ARCGIS !!!
L = (314376.57, 277097.9663, 253756.9869 ,265786.5632, 316712.6028, 232857.1468, 112063.9914, 135762.94, 131152.8206, 132323.662, 130317.865, 122075.3889, 19631.25154, 138142.7709, 199325.3382, 84158.4, 235101.9755, 179347.7901,
240720.646, 203442.04, 180101.061, 192130.6375, 243056.677, 159201.2212, 217305.5288, 62107.01542, 57496.89496, 58667.73631, 56661.93936, 48419.46322, 79510.13244, 99019.39624, 154407.9883, 95586.40056, 190184.6256, 134430.4403,
212852.124, 175573.5184, 152232.5389, 164262.1152, 215188.1548, 131332.6989, 189437.0064, 34238.49309, 29628.37263, 30799.21398, 28793.41703, 20761.78967, 115739.8845, 128061.9941, 151819.9124, 131816.1526, 187596.5497, 131842.3643,
156956.6541, 119678.0484, 97801.96823, 108366.6453, 159292.6849, 75437.22891, 131189.4431, 26207.94379, 30046.38946, 30022.33382, 36499.1228, 39154.16025, 165461.956, 173516.5609, 197274.4792, 181538.2241, 233051.1165, 177296.9311,
49979.50696, 36283.5152, 62559.4986, 17708.03221, 53338.17961, 39298.0475, 165299.2185, 133051.7993, 136890.245, 136866.1894, 143342.9784, 145998.0158, 272305.8115, 280360.4165, 304118.3347, 288382.0797, 339894.972, 284140.7867,
52935.33231, 128934.6503, 143847.8851, 99104.57653, 28403.13658, 116496.7582, 153039.7581, 210249.5101, 214087.9557, 214063.9001, 220540.6891, 223195.7265, 262646.9939, 357558.1272, 381316.0454, 365579.7904, 417092.6828,361338.4974,
342896.306, 305617.7003, 291163.1238, 294306.2972, 345232.3368, 261376.8808, 319481.1884, 164282.6751, 159672.5546, 160843.396, 158837.599, 150595.1229, 196032.549, 53287.46567, 18072.73241, 111682.1505, 47867.52331, 13615.3853,
205103.065, 190406.5166, 197459.3, 17835.09, 170439.194, 137982.189, 15022.6557, 168504.7907, 171893.2364, 1711869.181, 178345.9697, 181001.0072, 102411.6421, 206745.8721, 267147.8785, 174880.4788, 302924.5158, 247170.3304,
14088.98864, 175008.1185, 183668.6115, 162951.6989, 106219.1172, 123766.1215, 47423.29471, 188503.6722, 192342.1179, 192318.06, 198794.8512, 201449.8887, 157030.5306, 335812.2893, 359570.2076, 229499.3673, 395346.8449, 33959.65954,
193562.595, 156698.3255, 120368.67, 145386.92, 197099.806, 112302.1299, 174078.01, 22972.15, 22834.548, 22898.9, 20860.1683, 27742.407, 154363.5577, 139786.9169, 163892.2097, 179394.3768, 199558.9795, 144091.4906)
h= (75, 75, 75, 75, 75, 75, 1320, 75, 75, 75, 75, 75, 898, 75, 75, 665, 75, 75,
1369, 602, 349, 602, 1966, 778, 3837, 0 ,0, 0, 0, 0, 3384, 1602, 423, 2996, 423, 423,
1028, 260, 8, 260, 1624, 437, 3496, 0, 0, 0, 0, 0, 3043, 245, 0, 2655, 0, 0,
1822, 1054, 802, 1054, 2418, 1231, 4290, 0, 0, 0, 27, 0, 3837, 1039, 637, 3449, 637, 637,
893, 0, 0, 0, 1490, 0, 3362, 126, 126, 126, 126, 126, 2908, 126, 126, 2521, 126, 126,
161, 206, 306, 306, 306, 306, 1706, 306, 306, 306, 306, 306, 1325, 306, 306, 781, 306, 306,
1952, 1185, 287, 1185, 2549, 1361, 4421, 767, 767, 767, 767, 767, 3967, 1648, 137, 3674, 585, 0,
303, 758, 758, 758, 303, 758, 842, 758, 758, 758, 758, 758, 564, 168, 168, 1550, 168, 168,
2, 384, 384, 384, 2, 384, 1205, 384, 384, 384, 384, 384, 824, 384, 384, 1829, 384, 384,
1954, 1185, 746, 1185, 2520, 1362, 4401, 466, 466, 466, 437, 437, 3974, 1175, 319, 3790, 574, 319)
K = np.zeros(m*n)
N = np.zeros(m*n)
# In[61]:
#Definition of Objetive function
def objective(x):
global d
sum_obj = 0.0
for i in range(len(L)):
if x[i] < 0.01:
x[i]=0
if x[i] >= 0.276:
d= 1.484
elif 0.276 > x[i] >= 0.212:
d= 1.299
elif 0.212 > x[i] >= 0.148:
d= 1.086
elif 0.148 > x[i] >= 0.0975:
d=0.881
elif 0.0975 > x[i] >= 0.079:
d = 0.793
elif 0.079> x[i] >= 0.062:
d = 0.705
elif 0.062 > x[i] >= 0.049:
d= 0.626
elif 0.049 > x[i] >= 0.038:
d= 0.555
elif 0.038 > x[i] >= 0.030:
d=0.494
elif 0.030 > x[i] >= 0.024:
d=0.441
elif 0.024 > x[i] >= 0.0198:
d= 0.397
elif 0.0198 > x[i] >= 0.01565:
d= 0.353
elif 0.01565 > x[i] >= 0.0123:
d= 0.313
elif 0.0123 > x[i] >=0.0097:
d= 0.278
elif 0.0097 > x[i] >= 0.0077:
d=0.247
elif 0.0077 > x[i] >= 0.0061:
d=0.22
elif 0.0061 > x[i] >= 0.0049:
d=0.198
elif 0.0049 > x[i] >= 0.00389:
d=0.176
elif 0.00389 > x[i] >= 0.0032:
d= 0.159
elif 0.0032 > x[i] >= 0.0025 :
d= 0.141
elif 0.0025 > x[i] >= 0:
d= 0.123
K[i] = 10.674*L[i]/(C**1.852*d**4.871)
N[i] = K[i]*x[i]**1.852*J+x[i]*J*h[i]
sum_obj = sum_obj + N[i]
return sum_obj
#Definition of the constraints
GB=[]
for i in range(m):
GB.append(n*i*[0]+n*[1]+(m-1-i)*n*[0])
P=[]
for i in range(n):
P.append(m*(i*[0]+[1]+(n-1-i)*[0]))
DU=np.array(GB+P)
lb = np.array(m*[0] + demand) # Supply
ub = np.array(capacity + n*[np.inf]) # Demand
# In[62]:
# show initial objective
print('Initial SSE Objective: ' + str(objective(x0)))
# In[63]:
# optimize
bnds = Bounds(0.00, np.inf, keep_feasible=True)
cons = LinearConstraint(DU,lb,ub)
solution = minimize(objective,x0,method='SLSQP',constraints=cons,bounds=bnds,options={'iprint':2,'disp': True})
x = solution.x
Я не знаю почему, но в итоге моя модель дает меньше воды, чем нужно узлу. Какие-либо предложения ?