質問編集履歴

制約式が長いので大部分を省略しました。

2021/05/18 08:33

投稿

simpkins

スコア5

title CHANGED Viewed

File without changes

body CHANGED Viewed

@@ -135,36 +135,135 @@
 R_dem_OA = 0.2
 #関数
+def monte():
+    x = random.randrange(31, 61, 5)
+    Bcap_max = random.randrange(26, 176, 25)
+    Pb_max = random.randrange(17, 167, 25)
+    df=csv.resample('H').mean()
+    df = df.reset_index()
+    df['weekday'] =df['時間'].dt.weekday
+    df = df.set_index("時間")
+    df.loc[df['weekday'] < 5 , 'weekday'] = 1
+    df.loc[df['weekday'] >= 5, 'weekday'] = 0
+    def shuffle_days(df_month):
+        groups = [df for _, df in df_month.groupby('D')]
+        random.shuffle(groups)
+        return pd.concat(groups).reset_index(drop=True)
-    for i in range(1, T):
+    df = df.reset_index()
+    df['Y'] = df['時間'].dt.year
+    df['M'] = df['時間'].dt.month
+    df['D'] = df['時間'].dt.day
+    df = df.groupby(['Y', 'M','weekday']).apply(shuffle_days)
+    df = df.drop(['Y', 'M', 'D','weekday'], axis=1).reset_index(drop=True)
+    demand = list(df['demand'])
-        prb += Bcap[i] == Bcap[i-1] +  eff*Pb[i]
+    Ppv5 = list(df['pv_5'])
+    Cbuy = 15
+    d_term = [1,2,3]
+    tt = random.choice(d_term)
+    TT = 24 * tt
+    BBcap = [LpVariable('BBcap_{}'.format(i), None, None) for i in range(TT)]
+    PPb = [LpVariable('PPb_{}'.format(i), None, None) for i in range(TT)]
+    ddemand = [LpVariable('ddemand_{}'.format(i), 0, None) for i in range(TT)]
+    Psf = [LpVariable('Psf_{}'.format(i), None, None) for i in range(TT)]
+    Cbd = [LpVariable('Cbd_{}'.format(i), None, None) for i in range(TT)]
+    Csf = [LpVariable('Csf_{}'.format(i), None, None) for i in range(TT)]
+    Emp = [LpVariable('Emp_{}'.format(i), None, None) for i in range(TT)]
+    Eva = 10 #避難者数(Evacuee)
+）
+    x0 = [LpVariable('x0_{}'.format(i), 0, None) for i in range(TT)]
+    x1 = [LpVariable('x1_{}'.format(i), 0, None) for i in range(TT)]
+    x2 = [LpVariable('x2_{}'.format(i), 0, None) for i in range(TT)]
+    x3 = [LpVariable('x3_{}'.format(i), 0, None) for i in range(TT)]
+    SoC0_before = [0,0,0,0,0,0,0,0,0,0,0,0,0,0.3551,16.572,0,0,0,0,0,0,0,0,0]
+    SoC0_after = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,47.3011,50.3314,49.4792,0,0,0,0,0,0,0]
+    Rev_year = 4000000000
+    Cbuy = 15
+    Cp_end = LpVariable('Cp_end',None,None)
+    C0_before = 0.1
+    C0_after = 0
+    C1 = 0.5
+    C2 = 1
+    C3_before = 2
+    C3_after = 2.1
+    X0_before = 0.3
+    X0_after = 0.6 #コロナ禍
+    X1_before = 0.4
+    X1_after = 0.2 #コロナ禍
+    X2_before = 0.2
+    X2_after = 0.1 #コロナ禍
+    X3 = 0.1
+    rate_sf = 7
+    R_dem_OA = 0.2
+    date = random.randint(1, 365)
+    s = date * 24 + blackout_start_time
+    ddf2 = ddf[s:s+TT]
+    demand_usual_before = list(ddf2['demand_usual'])
+    demand_usual_after = list(ddf2['demand_usual_covid_19'])
+    PPpv5 = list(ddf2['pv_5'])
+    PPpv = [n * x for n in PPpv5] #PV発電量
+    demand_hour_Emp = list(ddf2['demand_hour_Emp'])
+    Emp_usual_before = list(ddf2['Emp_usual'])
+    Emp_usual_after = list(ddf2['Emp_usual_covid_19'])
+    PPpv = [n * x for n in PPpv5]
+    BL = list(ddf2['Base'])
+    Emp_total = sum(ddf['Emp_usual'])
+    Rev_hour_Emp = Rev_year / Emp_total #
+    if covid_19 == 0
+        demand_usual = demand_usual_before
+        Emp_usual = Emp_usual_before
+        SoC0 = SoC0_before[blackout_start_time] / 100
+        C0 = C0_before
+        C3 = C3_before
+        X0 = X0_before
+        X1 = X1_before
+        X2 = X2_before
+    else:
+        demand_usual = demand_usual_after
+        Emp_usual = Emp_usual_after
+        SoC0 = SoC0_after[blackout_start_time] / 100
+        C0 = C0_after
+        C3 = C3_after
+        X0 = X0_after
+        X1 = X1_after
+        X2 = X2_after
+    CX_elec = C0 * X0 + C1 * X1 + C2 * X2 + C3 * X3
+# Optimisation problem　目的関数(最小化)
+    prb = LpProblem('Battery Operation')
+# Objective　目的関数
+    objective1= lpSum(Cp) +Binterest*(CCB_inv*Bcap_max+CPB_inv*Pb_max)+PVinterest*Cpv_inv*Ppv_rated +1700*12*Pgrid_max
+    objective2 = Cp_end + lpSum(Cbd) + lpSum(Csf)
+    prb += objective1 + objective2
+    # Constraints　制約
+    prb += Ppv_rated == 9*x
+    for i in range(T):
+        if Pgrid[i]>=0:
+            Cp[i] == Cbuy * Pgrid[i]
+        else:
+            Cp[i] == Csell * Pgrid[i]
+        prb+= Ppv[i]== Ppv5[i]* x
 ###ここからは災害の制約
     for i in range(TT):
         prb += Cbd[i] == Rev_hour_Emp * Emp_usual_before[i] * (CX_elec - (C0 * x0[i])\
         - (C1 * x1[i]) - (C2 * x2[i]) - (C3 * x3[i]))
-        prb += Csf[i] == rate_sf * Rev_hour_Emp * Psf[i] / demand_hour_Emp[i]
-        prb += 0 <= x0[i] <= X0
-        prb += 0 <= x1[i] <= X1
-        prb += 0 <= x2[i] <= X2
-        prb += 0 <= x3[i] <= X3
-        prb += Emp[i] == Emp_usual_before[i] * (x0[i] + x1[i] + x2[i] + x3[i])
-        prb += ddemand[i] == BL[i] + demand_hour_Emp[i] * (Emp[i] + (1 - R_dem_OA) * Eva)
-        prb += BBcap[i] <= Bcap_max
-        prb += BBcap[i] >= 0
-        prb += PPb[i] <= Pb_max
-        prb += PPb[i] >= (-1) * Pb_max
-        prb += Psf[i] >= 0
-        prb += PPb[i] <= PPpv[i] + Psf[i] - ddemand[i]
-        if i == 0:
-            prb += BBcap[0] == SoC0 * Bcap_max + eff * PPb[0]
-        else:
-            prb += BBcap[i] == BBcap[i-1] + eff * PPb[i]
-        prb += Cp_end == Cbuy * (SoC0 * Bcap_max - BBcap[TT-1])
     prb.solve()
     cost = pulp.value(objective1)