# Given price function
def price_function(X, Y, K_upper, K_lower, theta, B):
ratio = Y/X
return K_upper - (K_upper - K_lower) / (1 + np.exp(-B * (ratio - theta)))
# Objective function for joint optimization of B and theta
def optimization_objective(params, Y_val, K_upper_val, K_lower_val, A):
B_val, theta_val = params
# The price when demand is exactly at the lower transition (K_lower * Y_val)
price_at_lower_transition = price_function(Y_val * K_lower_val, Y_val, K_upper_val, K_lower_val, theta_val, B_val)
# The price when demand is exactly at the upper transition (A * K_upper_val)
price_at_upper_transition = price_function(A * Y_val * K_upper_val, Y_val, K_upper_val, K_lower_val, theta_val, B_val)
# Ensure the price at the lower transition is close to K_lower and at the upper transition is close to K_upper
term1 = (price_at_lower_transition - K_lower_val)**2
term2 = (price_at_upper_transition - K_upper_val)**2
return term1 + term2
def market_welfare_for_A(A, consumption, production, K_upper, K_lower):
total_producer_welfare_per_unit = 0
total_consumer_welfare_per_unit = 0
total_producer_welfare = 0
total_consumer_welfare = 0
total_units_transacted = 0
for c, p in zip(consumption, production):
# Run the optimization for B and theta given A
initial_guess = [2, 0.5] # [B, theta]
bnds = [(0, 100), (0, 100)]
# Minimize
result = minimize(optimization_objective, initial_guess, args=(1, K_upper_val, K_lower_val, A), method='L-BFGS-B', bounds=bnds)
optimal_B, optimal_theta = result.x
# Price at this timestep for the given A
price = price_function(c, p, K_upper, K_lower, optimal_theta, optimal_B)
# because there demand is always much larger
units_transacted = p
total_units_transacted += units_transacted
# Producer welfare per unit: revenue from sold energy
producer_welfare_per_unit =(units_transacted * price - units_transacted * K_lower) if units_transacted != 0 else 0
# Consumer welfare per unit: total value to consumers minus what they paid, per unit
consumer_welfare_per_unit = (units_transacted * K_upper - units_transacted * price) if units_transacted != 0 else 0
total_producer_welfare += producer_welfare_per_unit
total_consumer_welfare += consumer_welfare_per_unit
total_units_transacted += units_transacted
# Average welfare per unit transacted for the day (producer and consumer)
average_producer_welfare_per_unit = total_producer_welfare / total_units_transacted
average_consumer_welfare_per_unit = total_consumer_welfare / total_units_transacted
# Min-max welfare metric
min_welfare_per_unit = min(average_producer_welfare_per_unit, average_consumer_welfare_per_unit)
return -min_welfare_per_unit, average_producer_welfare_per_unit, average_consumer_welfare_per_unit
def objective_for_A(A, consumption, production, K_upper, K_lower):
welfare_value, _, _ = market_welfare_for_A(A[0], consumption, production, K_upper, K_lower)
return welfare_value
# Example data
consumption = demands # for example
production = supplies # for example
K_upper_val = 0.34
K_lower_val = 0.08
# Bounds for A
bounds = [(0, 500)]
# Use differential_evolution to optimize A
result = differential_evolution(objective_for_A, bounds, args=(consumption, production, K_upper_val, K_lower_val))
optimal_A = result.x[0]
best_welfare_value, best_producer_welfare, best_consumer_welfare = market_welfare_for_A(optimal_A, consumption, production, K_upper_val, K_lower_val)
print(f"Optimal A: {optimal_A}")
print(f"Best average welfare: {best_welfare_value}")
print(f"Best producer welfare: {best_producer_welfare}")
print(f"Best consumer welfare: {best_consumer_welfare}")