File: test3.py

package info (click to toggle)
python-pulp 2.6.0%2Bdfsg-1
  • links: PTS, VCS
  • area: main
  • in suites: bookworm
  • size: 14,720 kB
  • sloc: python: 7,505; makefile: 16; sh: 16
file content (131 lines) | stat: -rw-r--r-- 3,732 bytes parent folder | download | duplicates (2)
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
#!/usr/bin/env python
# @(#) $Jeannot: test3.py,v 1.3 2004/03/20 17:06:54 js Exp $

# Deterministic generation planning using mixed integer linear programming.

# The goal is to minimise the cost of generation while satisfaying demand
# using a few thermal units and an hydro unit.
# The thermal units have a proportional cost and a startup cost.
# The hydro unit has an initial storage.

from __future__ import print_function
from pulp import *
from math import *

prob = LpProblem("test3", LpMinimize)

# The number of time steps
tmax = 9
# The number of thermal units
units = 5
# The minimum demand
dmin = 10.0
# The maximum demand
dmax = 150.0
# The maximum thermal production
tpmax = 150.0
# The maximum hydro production
hpmax = 100.0
# Initial hydro storage
sini = 50.0

# Time range
time = list(range(tmax))
# Time range (and one more step for the last state of plants)
xtime = list(range(tmax + 1))
# Units range
unit = list(range(units))
# The demand
demand = [
    dmin + (dmax - dmin) * 0.5 + 0.5 * (dmax - dmin) * sin(4 * t * 2 * 3.1415 / tmax)
    for t in time
]
# Maximum output for the thermal units
pmax = [tpmax / units for i in unit]
# Minimum output for the thermal units
pmin = [tpmax / (units * 3.0) for i in unit]
# Proportional cost of the thermal units
costs = [i + 1 for i in unit]
# Startup cost of the thermal units.
startupcosts = [100 * (i + 1) for i in unit]

# Production variables for each time step and each thermal unit.
p = LpVariable.matrix("p", (time, unit), 0)
for t in time:
    for i in unit:
        p[t][i].upBound = pmax[i]

# State (started/stopped) variables for each time step and each thermal unit
d = LpVariable.matrix("d", (xtime, unit), 0, 1, LpInteger)

# Production constraint relative to the unit state (started/stoped)
for t in time:
    for i in unit:
        # If the unit is not started (d==0) then p<=0 else p<=pmax
        prob += p[t][i] <= pmax[i] * d[t][i]
        # If the unit is not started then p>=0 else p>= pmin
        prob += p[t][i] >= pmin[i] * d[t][i]

# Startup variables: 1 if the unit will be started next time step
u = LpVariable.matrix("u", (time, unit), 0)

# Dynamic startup constraints
# Initialy, all groups are started
for t in time:
    for i in unit:
        # u>=1 if the unit is started next time step
        prob += u[t][i] >= d[t + 1][i] - d[t][i]

# Storage for the hydro plant (must not go below 0)
s = LpVariable.matrix("s", xtime, 0)

# Initial storage
s[0] = sini

# Hydro production
ph = [s[t] - s[t + 1] for t in time]
for t in time:
    # Must be positive (no pumping)
    prob += ph[t] >= 0
    # And lower than hpmax
    prob += ph[t] <= hpmax

# Total production must equal demand
for t in time:
    prob += demand[t] == lpSum(p[t]) + ph[t]

# Thermal production cost
ctp = lpSum([lpSum([p[t][i] for t in time]) * costs[i] for i in unit])
# Startup costs
cts = lpSum([lpSum([u[t][i] for t in time]) * startupcosts[i] for i in unit])
# The objective is the total cost
prob += ctp + cts

# Solve the problem
prob.solve()

print("Minimum total cost:", prob.objective.value())

# Print the results
print("   D    S     U ", end=" ")
for i in unit:
    print("  T%d    " % i, end=" ")
print()

for t in time:
    # Demand, hydro storage, hydro production
    print("%5.1f" % demand[t], "%5.1f" % value(s[t]), "%5.1f" % value(ph[t]), end=" ")
    for i in unit:
        # Thermal production
        print("%4.1f" % value(p[t][i]), end=" ")
        # The state of the unit
        if value(d[t][i]):
            print("+", end=" ")
        else:
            print("-", end=" ")
        # Wether the unit will be started
        if value(u[t][i]):
            print("*", end=" ")
        else:
            print(" ", end=" ")
    print()