質問編集履歴

1

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

2021/05/18 08:33

投稿

simpkins
simpkins

スコア5

test CHANGED
File without changes
test CHANGED
@@ -272,13 +272,251 @@
272
272
 
273
273
  #関数
274
274
 
275
-
275
+ def monte():
276
+
276
-
277
+ x = random.randrange(31, 61, 5)
278
+
279
+ Bcap_max = random.randrange(26, 176, 25)
280
+
281
+ Pb_max = random.randrange(17, 167, 25)
282
+
283
+
284
+
285
+ df=csv.resample('H').mean()
286
+
287
+ df = df.reset_index()
288
+
289
+ df['weekday'] =df['時間'].dt.weekday
290
+
291
+ df = df.set_index("時間")
292
+
293
+ df.loc[df['weekday'] < 5 , 'weekday'] = 1
294
+
295
+ df.loc[df['weekday'] >= 5, 'weekday'] = 0
296
+
297
+
298
+
299
+ def shuffle_days(df_month):
300
+
301
+ groups = [df for _, df in df_month.groupby('D')]
302
+
303
+ random.shuffle(groups)
304
+
305
+ return pd.concat(groups).reset_index(drop=True)
306
+
307
+
308
+
309
+ df = df.reset_index()
310
+
311
+ df['Y'] = df['時間'].dt.year
312
+
313
+ df['M'] = df['時間'].dt.month
314
+
315
+ df['D'] = df['時間'].dt.day
316
+
317
+ df = df.groupby(['Y', 'M','weekday']).apply(shuffle_days)
318
+
319
+ df = df.drop(['Y', 'M', 'D','weekday'], axis=1).reset_index(drop=True)
320
+
321
+
322
+
323
+ demand = list(df['demand'])
324
+
325
+ Ppv5 = list(df['pv_5'])
326
+
327
+ Cbuy = 15
328
+
329
+ d_term = [1,2,3]
330
+
331
+ tt = random.choice(d_term)
332
+
333
+
334
+
335
+ TT = 24 * tt
336
+
337
+ BBcap = [LpVariable('BBcap_{}'.format(i), None, None) for i in range(TT)]
338
+
339
+ PPb = [LpVariable('PPb_{}'.format(i), None, None) for i in range(TT)]
340
+
341
+ ddemand = [LpVariable('ddemand_{}'.format(i), 0, None) for i in range(TT)]
342
+
343
+ Psf = [LpVariable('Psf_{}'.format(i), None, None) for i in range(TT)]
344
+
345
+
346
+
347
+ Cbd = [LpVariable('Cbd_{}'.format(i), None, None) for i in range(TT)]
348
+
349
+ Csf = [LpVariable('Csf_{}'.format(i), None, None) for i in range(TT)]
350
+
351
+
352
+
353
+ Emp = [LpVariable('Emp_{}'.format(i), None, None) for i in range(TT)]
354
+
355
+ Eva = 10 #避難者数(Evacuee)
356
+
357
+
358
+
359
+ x0 = [LpVariable('x0_{}'.format(i), 0, None) for i in range(TT)]
360
+
361
+ x1 = [LpVariable('x1_{}'.format(i), 0, None) for i in range(TT)]
362
+
363
+ x2 = [LpVariable('x2_{}'.format(i), 0, None) for i in range(TT)]
364
+
365
+ x3 = [LpVariable('x3_{}'.format(i), 0, None) for i in range(TT)]
366
+
367
+
368
+
369
+ 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]
370
+
371
+ 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]
372
+
373
+
374
+
375
+ Rev_year = 4000000000
376
+
377
+ Cbuy = 15
378
+
379
+ Cp_end = LpVariable('Cp_end',None,None)
380
+
381
+ C0_before = 0.1
382
+
383
+ C0_after = 0
384
+
385
+ C1 = 0.5
386
+
387
+ C2 = 1
388
+
389
+ C3_before = 2
390
+
391
+ C3_after = 2.1
392
+
393
+
394
+
395
+ X0_before = 0.3
396
+
397
+ X0_after = 0.6 #コロナ禍
398
+
399
+ X1_before = 0.4
400
+
401
+ X1_after = 0.2 #コロナ禍
402
+
403
+ X2_before = 0.2
404
+
405
+ X2_after = 0.1 #コロナ禍
406
+
407
+ X3 = 0.1
408
+
409
+ rate_sf = 7
410
+
411
+ R_dem_OA = 0.2
412
+
413
+ date = random.randint(1, 365)
414
+
415
+ s = date * 24 + blackout_start_time
416
+
417
+ ddf2 = ddf[s:s+TT]
418
+
419
+
420
+
421
+ demand_usual_before = list(ddf2['demand_usual'])
422
+
423
+ demand_usual_after = list(ddf2['demand_usual_covid_19'])
424
+
425
+
426
+
427
+ PPpv5 = list(ddf2['pv_5'])
428
+
429
+ PPpv = [n * x for n in PPpv5] #PV発電量
430
+
431
+
432
+
433
+ demand_hour_Emp = list(ddf2['demand_hour_Emp'])
434
+
435
+ Emp_usual_before = list(ddf2['Emp_usual'])
436
+
437
+ Emp_usual_after = list(ddf2['Emp_usual_covid_19'])
438
+
439
+ PPpv = [n * x for n in PPpv5]
440
+
441
+
442
+
443
+ BL = list(ddf2['Base'])
444
+
445
+ Emp_total = sum(ddf['Emp_usual'])
446
+
447
+ Rev_hour_Emp = Rev_year / Emp_total #
448
+
449
+
450
+
451
+ if covid_19 == 0
452
+
453
+ demand_usual = demand_usual_before
454
+
455
+ Emp_usual = Emp_usual_before
456
+
457
+ SoC0 = SoC0_before[blackout_start_time] / 100
458
+
459
+ C0 = C0_before
460
+
461
+ C3 = C3_before
462
+
463
+ X0 = X0_before
464
+
465
+ X1 = X1_before
466
+
467
+ X2 = X2_before
468
+
469
+ else:
470
+
471
+ demand_usual = demand_usual_after
472
+
473
+ Emp_usual = Emp_usual_after
474
+
475
+ SoC0 = SoC0_after[blackout_start_time] / 100
476
+
477
+ C0 = C0_after
478
+
479
+ C3 = C3_after
480
+
481
+ X0 = X0_after
482
+
483
+ X1 = X1_after
484
+
485
+ X2 = X2_after
486
+
487
+
488
+
489
+ CX_elec = C0 * X0 + C1 * X1 + C2 * X2 + C3 * X3
490
+
491
+ # Optimisation problem 目的関数(最小化)
492
+
493
+ prb = LpProblem('Battery Operation')
494
+
495
+ # Objective 目的関数
496
+
497
+ objective1= lpSum(Cp) +Binterest*(CCB_inv*Bcap_max+CPB_inv*Pb_max)+PVinterest*Cpv_inv*Ppv_rated +1700*12*Pgrid_max
498
+
499
+ objective2 = Cp_end + lpSum(Cbd) + lpSum(Csf)
500
+
501
+ prb += objective1 + objective2
502
+
503
+ # Constraints 制約
504
+
505
+ prb += Ppv_rated == 9*x
506
+
277
- for i in range(1, T):
507
+ for i in range(T):
508
+
278
-
509
+ if Pgrid[i]>=0:
510
+
511
+ Cp[i] == Cbuy * Pgrid[i]
512
+
513
+ else:
514
+
515
+ Cp[i] == Csell * Pgrid[i]
516
+
517
+
518
+
279
- prb += Bcap[i] == Bcap[i-1] + eff*Pb[i]
519
+ prb+= Ppv[i]== Ppv5[i]* x
280
-
281
-
282
520
 
283
521
  ###ここからは災害の制約
284
522
 
@@ -290,47 +528,7 @@
290
528
 
291
529
  - (C1 * x1[i]) - (C2 * x2[i]) - (C3 * x3[i]))
292
530
 
293
-
294
-
295
- prb += Csf[i] == rate_sf * Rev_hour_Emp * Psf[i] / demand_hour_Emp[i]
531
+
296
-
297
- prb += 0 <= x0[i] <= X0
298
-
299
- prb += 0 <= x1[i] <= X1
300
-
301
- prb += 0 <= x2[i] <= X2
302
-
303
- prb += 0 <= x3[i] <= X3
304
-
305
- prb += Emp[i] == Emp_usual_before[i] * (x0[i] + x1[i] + x2[i] + x3[i])
306
-
307
-
308
-
309
- prb += ddemand[i] == BL[i] + demand_hour_Emp[i] * (Emp[i] + (1 - R_dem_OA) * Eva)
310
-
311
- prb += BBcap[i] <= Bcap_max
312
-
313
- prb += BBcap[i] >= 0
314
-
315
- prb += PPb[i] <= Pb_max
316
-
317
- prb += PPb[i] >= (-1) * Pb_max
318
-
319
- prb += Psf[i] >= 0
320
-
321
- prb += PPb[i] <= PPpv[i] + Psf[i] - ddemand[i]
322
-
323
- if i == 0:
324
-
325
- prb += BBcap[0] == SoC0 * Bcap_max + eff * PPb[0]
326
-
327
- else:
328
-
329
- prb += BBcap[i] == BBcap[i-1] + eff * PPb[i]
330
-
331
-
332
-
333
- prb += Cp_end == Cbuy * (SoC0 * Bcap_max - BBcap[TT-1])
334
532
 
335
533
  prb.solve()
336
534