Define the power system parameters

num_generators = 3 num_buses = 4 power_demand = [200, 150, 100, 250] power_generation_min = [50, 40, 30] power_generation_max = [150, 120, 100] generation_cost = [10, 12, 15] bus_voltage_min = [0.95, 0.95, 0.95, 0.95] bus_voltage_max = [1.05, 1.05, 1.05, 1.05] line_flow_max = [100, 100, 100, 100]

Create the optimization model

model = ConcreteModel()

Define the decision variables

model.power_generation = Var(range(num_generators), within=NonNegativeReals)

Define the objective function

model.cost = Objective(expr=sum(generation_cost[i] * model.power_generation[i] for i in range(num_generators)))

Define the power balance constraint

def power_balance_rule(model, i): return sum(model.power_generation[j] for j in range(num_generators)) == power_demand[i] model.power_balance = Constraint(range(num_buses), rule=power_balance_rule)

Define the power generation limits constraint

def power_generation_limits_rule(model, i): return power_generation_min[i] <= model.power_generation[i] <= power_generation_max[i] model.power_generation_limits = Constraint(range(num_generators), rule=power_generation_limits_rule)

Define the bus voltage limits constraint

def bus_voltage_limits_rule(model, i): return bus_voltage_min[i] <= sum(model.power_generation[j] for j in range(num_generators)) <= bus_voltage_max[i] model.bus_voltage_limits = Constraint(range(num_buses), rule=bus_voltage_limits_rule)

Define the line flow limits constraint

def line_flow_limits_rule(model, i): return sum(model.power_generation[j] for j in range(num_generators)) <= line_flow_max[i] model.line_flow_limits = Constraint(range(num_buses), rule=line_flow_limits_rule)

Solve the optimization problem

solver = SolverFactory('glpk') solver.solve(model)

Print the optimal solution

print("Optimal power generation:") for i in range(num_generators): print(f"Generator {i+1}: {model.power_generation[i].value}"