Skip to content
AVOA_UCP_NEW
import numpy as np
import random
import matplotlib.pyplot as plt
# ------------------------------
# 1. SYSTEM PARAMETERS
# ------------------------------
num_units = 10
num_hours = 24
# Initial ON/OFF status (hrs ON >0, OFF <0)
initial_status = np.array([ 8, 8, 8, -8, -8, -8, 8, -8, -8, -8])
# Generation limits
Pmin = np.array([150,150,20,20,25,20,25,10,10,10])
Pmax = np.array([455,455,130,130,162,80,85,55,55,55])
# Fuel cost coefficients
a = np.array([0.00048,0.00031,0.00200,0.00211,0.00398,0.00712,0.00079,0.00413,0.002202,0.00173])
b = np.array([16.19,17.26,16.60,16.50,19.70,22.26,27.74,25.92,27.27,27.79])
c = np.array([1000,970,700,680,450,370,480,660,665,670])
# Startup costs
SU_hot = np.array([4500,5000,550,560,900,170,260,30,30,30])
SU_cold = np.array([9000,10000,1100,1120,1800,340,520,60,60,60])
T_cold = np.array([5,5,4,4,4,2,2,0,0,0])
# Min up/down times
T_MU = np.array([8,8,5,5,6,3,3,1,1,1])
T_MD = np.array([8,8,5,5,6,3,3,1,1,1])
# Load demand
load = np.array([700,750,850,950,1000,1100,1150,1200,
1300,1400,1450,1500,1400,1300,1200,1050,
1000,1100,1200,1400,1300,1100,900,800])
reserve_pct = 0.10 # 10% spinning reserve
# AVOA parameters
pop_size = 30
max_iter = 300
# ------------------------------
# 2. GREEDY SEED
# ------------------------------
def greedy_seed():
"""Heuristic: for each hour, turn on cheapest units until load+reserve met."""
seed = np.zeros((num_units, num_hours), int)
cost_rank = a*(Pmax**2) + b*Pmax + c # approximate full-load cost
order = np.argsort(cost_rank)
for t in range(num_hours):
target = load[t]*(1+reserve_pct)
cum = 0
for i in order:
seed[i,t] = 1
cum += Pmin[i]
if cum >= target:
break
return seed
# ------------------------------
# 3. POPULATION INIT
# ------------------------------
def init_pop():
pop = [greedy_seed()]
for _ in range(pop_size-1):
pop.append((np.random.rand(num_units,num_hours)>0.5).astype(int))
return pop# ------------------------------
# 4. LAMBDA DISPATCH
# ------------------------------
def lambda_dispatch(schedule):
dispatch = np.zeros_like(schedule, float)
for t in range(num_hours):
eff_load = load[t]*(1+reserve_pct)
on = np.where(schedule[:,t]==1)[0]
if len(on)==0: continue
lam_lo, lam_hi = 0.0, 1e4
for _ in range(100):
lam = 0.5*(lam_lo+lam_hi)
gen = (lam - b[on])/(2*a[on])
gen = np.clip(gen, Pmin[on], Pmax[on])
if abs(gen.sum()-eff_load)<1e-2: break
if gen.sum()>eff_load: lam_hi=lam
else: lam_lo=lam
dispatch[on,t] = gen
dispatch_matrix = np.round(dispatch)
return dispatch, dispatch_matrix
# ------------------------------
# 5. FITNESS EVALUATION
# ------------------------------
def evaluate(schedule, pen_power, pen_up):
status = (initial_status>0).astype(int)
ton = np.where(initial_status>0, initial_status, 0).tolist()
toff = np.where(initial_status<0, -initial_status, 0).tolist()
cost = 0.0
disp = lambda_dispatch(schedule)
for t in range(num_hours):
eff = load[t]*(1+reserve_pct)
total = disp[:,t].sum()
if abs(total-eff)>1e-2:
cost += pen_power
for i in range(num_units):
u = schedule[i,t]
p = disp[i,t]
if u:
cost += a[i]*p*p + b[i]*p + c[i]
if status[i]==0:
cost += SU_hot[i] if toff[i]<=T_cold[i] else SU_cold[i]
ton[i] = ton[i]+1 if status[i] else 1
toff[i] = 0
else:
if status[i] and ton[i]<T_MU[i]:
cost += pen_up
toff[i] = toff[i]+1 if status[i]==0 else 1
ton[i] = 0
status[i]=u
return cost
# ------------------------------
# 6. LOCAL SEARCH (MEMETIC)
# ------------------------------
def local_search(sol, best_cost):
best = sol.copy()
bcost = best_cost
for _ in range(100):
i = random.randrange(num_units)
t = random.randrange(num_hours)
cand = best.copy()
cand[i,t] ^= 1
c = evaluate(cand, pen_power=1e6, pen_up=1e5)
if c<bcost:
best, bcost = cand, c
return best, bcost
# ------------------------------
# 7. EXPLORATION UPDATE
# ------------------------------
def update_loc(current, best, G):
r = np.random.rand(*current.shape)
C = np.abs(2*r + best - current)
new = best - G*C
return (new>0.5).astype(int)
# ------------------------------
# 8. AVOA MAIN
# ------------------------------
def avoa_ucp():
pop = init_pop()
# evaluate initial
fitness = [evaluate(sol,1e6,1e5) for sol in pop]
best_idx = np.argmin(fitness)
best_sol, best_cost = pop[best_idx], fitness[best_idx]
hist = []
G0 = 2.0
for it in range(max_iter):
G = G0*(1 - it/max_iter)
pen_p = 1e6*(1+it/max_iter)
pen_u = 1e5*(1+it/max_iter)
new_pop = []
for sol in pop:
# update by exploration
sol = update_loc(sol, best_sol, G)
new_pop.append(sol)
pop = new_pop
# re-evaluate
for sol in pop:
f = evaluate(sol, pen_p, pen_u)
if f<best_cost:
best_cost, best_sol = f, sol.copy()
# memetic polishing
if it%20==0:
best_sol, best_cost = local_search(best_sol, best_cost)
hist.append(best_cost)
if it%20==0:
print(f"Iter {it}: Cost = {best_cost:.2f}")
# plot
plt.plot(hist)
plt.xlabel('Iter')
plt.ylabel('Best Cost')
plt.title('AVOA UCP')
plt.grid(True)
plt.show()
return best_sol, best_cost
# ------------------------------
# 9. RUN
# ------------------------------
best_schedule, best_cost = avoa_ucp()
print("\nFinal Best Cost:", best_cost)
print("Schedule matrix (unit × hour):\n", best_schedule)dispatch, dispatch_matrix = lambda_dispatch(best_schedule)
np.transpose(dispatch_matrix)def calculate_total_cost(schedule, dispatch):
total_cost = 0
status_prev = [1 if h > 0 else 0 for h in initial_status]
T_on = [h if h > 0 else 0 for h in initial_status]
T_off = [-h if h < 0 else 0 for h in initial_status]
for t in range(num_hours):
for i in range(num_units):
u = schedule[t][i]
p = dispatch[t][i]
if u == 1:
# Fuel cost
total_cost += a[i] * (p ** 2) + b[i] * p + c[i]
# Startup cost
if status_prev[i] == 0:
if T_off[i] <= T_cold[i]:
total_cost += SU_hot[i]
else:
total_cost += SU_cold[i]
# Update ON/OFF duration
T_on[i] = T_on[i] + 1 if status_prev[i] == 1 else 1
T_off[i] = 0
else:
if status_prev[i] == 1 and T_on[i] < T_MU[i]:
total_cost += 1e5 # Penalty for violating min up time
T_off[i] = T_off[i] + 1 if status_prev[i] == 0 else 1
T_on[i] = 0
status_prev[i] = u
return total_costtotal_cost = calculate_total_cost(np.transpose(best_schedule),np.transpose(dispatch_matrix))
total_cost