Untitled

# 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}")