Skip to content

Commit 75c43f8

Browse files
TransurgeonTransurgeon
committed
add examples for power systems
1 parent 263117b commit 75c43f8

File tree

6 files changed

+269
-2
lines changed

6 files changed

+269
-2
lines changed

‎cvxpy/sandbox/bilinear.py‎

Lines changed: 41 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,41 @@
1+
#!/usr/bin/env python3
2+
# CVXPY implementation of the bilinear optimization problem
3+
# This example assumes CVXPY can handle non-convex bilinear constraints
4+
5+
import cvxpy as cp
6+
7+
# First version: continuous variables
8+
print("=" * 50)
9+
print("Version 1: Continuous variables")
10+
print("=" * 50)
11+
12+
# Create variables (non-negative)
13+
x = cp.Variable(nonneg=True, name="x")
14+
y = cp.Variable(nonneg=True, name="y")
15+
z = cp.Variable(nonneg=True, name="z")
16+
17+
# Set objective: maximize x
18+
objective = cp.Maximize(x)
19+
20+
# Define constraints
21+
constraints = [
22+
# Linear constraint: x + y + z <= 10
23+
x + y + z <= 10,
24+
25+
# Bilinear inequality constraint: x * y <= 2
26+
x * y <= 2,
27+
28+
# Bilinear equality constraint: x * z + y * z == 1
29+
x * z + y * z == 1
30+
]
31+
32+
# Create and solve problem
33+
problem = cp.Problem(objective, constraints)
34+
problem.solve(solver=cp.IPOPT, verbose=True, nlp=True)
35+
36+
# Print results
37+
print(f"Status: {problem.status}")
38+
print(f"Optimal value: {problem.value:.6f}")
39+
print(f"x = {x.value:.6f}")
40+
print(f"y = {y.value:.6f}")
41+
print(f"z = {z.value:.6f}")

‎cvxpy/sandbox/bilinear_gurobi.py‎

Lines changed: 44 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,44 @@
1+
#!/usr/bin/env python3.11
2+
3+
# Copyright 2025, Gurobi Optimization, LLC
4+
5+
# This example formulates and solves the following simple bilinear model:
6+
# maximize x
7+
# subject to x + y + z <= 10
8+
# x * y <= 2 (bilinear inequality)
9+
# x * z + y * z = 1 (bilinear equality)
10+
# x, y, z non-negative (x integral in second version)
11+
12+
import gurobipy as gp
13+
from gurobipy import GRB
14+
15+
# Create a new model
16+
m = gp.Model("bilinear")
17+
18+
# Create variables
19+
x = m.addVar(name="x")
20+
y = m.addVar(name="y")
21+
z = m.addVar(name="z")
22+
23+
# Set objective: maximize x
24+
m.setObjective(1.0 * x, GRB.MAXIMIZE)
25+
26+
# Add linear constraint: x + y + z <= 10
27+
m.addConstr(x + y + z <= 10, "c0")
28+
29+
# Add bilinear inequality constraint: x * y <= 2
30+
m.addConstr(x * y <= 2, "bilinear0")
31+
32+
# Add bilinear equality constraint: x * z + y * z == 1
33+
m.addConstr(x * z + y * z == 1, "bilinear1")
34+
35+
# Optimize model
36+
m.optimize()
37+
38+
m.printAttr("x")
39+
40+
# Constrain 'x' to be integral and solve again
41+
x.VType = GRB.INTEGER
42+
m.optimize()
43+
44+
m.printAttr("x")

‎cvxpy/sandbox/control_of_car_matrix.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -4,7 +4,7 @@
44
import cvxpy as cp
55

66

7-
def solve_car_control(x_final, L=0.1, N=2, h=0.1, gamma=10):
7+
def solve_car_control(x_final, L=0.1, N=30, h=0.1, gamma=10):
88
"""
99
Solve the nonlinear optimal control problem for car trajectory planning.
1010

‎cvxpy/sandbox/gurobipy_acopf.py‎

Lines changed: 107 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,107 @@
1+
#!/usr/bin/env python
2+
3+
# Copyright 2025, Gurobi Optimization, LLC
4+
5+
# A Power Flow Problem is illustrated on a 4 bus (node) network using an AC (Alternating Current) represention
6+
#
7+
#
8+
# Zij = Impedance between i and j P_i = Active Power Injected into bus i
9+
# Yij = Admittance between i and j (Yij = 1/Zij) Q_i = Reactive Power Injeced into bus i
10+
# Yij = Gij + j Bij (j = sqrt(-1)) V_i = Voltage Magnitude at bus i
11+
# Theta_i = Voltage Angle at bus i
12+
#
13+
# Real (P) and Reactive (Q) Power balance equations must be met on all buses
14+
#
15+
# P_i = V_i \sum_{j=1}^{N} V_j (G_{ij} \cos(\theta_i - \theta_j) + B_{ij} \sin(\theta_i - \theta_j))
16+
#
17+
# Q_i = V_i \sum_{j=1}^{N} V_j (G_{ij} \sin(\theta_i - \theta_j) - B_{ij} \cos(\theta_i - \theta_j))
18+
#
19+
#
20+
# (1)---(2) Bus 1: Swing Bus. Fixed Voltage Magnitude and Angle. Variable P and Q.
21+
# | | Bus 2: Load Bus. Fixed P and Q. Variable Voltage Magnitude and Angle.
22+
# | | Bus 3. Generation bus. Fixed Voltage Magnitude and P. Variable Voltage Angle and Q.
23+
# (4)---(3) Bus 4. Load Bus. Fixed P and Q. Variable Voltage and Angle.
24+
#
25+
# Objective: minimize overall reactive power on buses 1 and 3
26+
import numpy as np
27+
import gurobipy as gp
28+
from gurobipy import GRB
29+
from gurobipy import nlfunc
30+
31+
32+
# Number of Buses (Nodes)
33+
N = 4
34+
35+
# Conductance/susceptance components
36+
G = np.array(
37+
[
38+
[1.7647, -0.5882, 0.0, -1.1765],
39+
[-0.5882, 1.5611, -0.3846, -0.5882],
40+
[0.0, -0.3846, 1.5611, -1.1765],
41+
[-1.1765, -0.5882, -1.1765, 2.9412],
42+
]
43+
)
44+
B = np.array(
45+
[
46+
[-7.0588, 2.3529, 0.0, 4.7059],
47+
[2.3529, -6.629, 1.9231, 2.3529],
48+
[0.0, 1.9231, -6.629, 4.7059],
49+
[4.7059, 2.3529, 4.7059, -11.7647],
50+
]
51+
)
52+
53+
# Assign bounds where fixings are needed
54+
v_lb = np.array([1.0, 0.0, 1.0, 0.0])
55+
v_ub = np.array([1.0, 1.5, 1.0, 1.5])
56+
P_lb = np.array([-3.0, -0.3, 0.3, -0.2])
57+
P_ub = np.array([3.0, -0.3, 0.3, -0.2])
58+
Q_lb = np.array([-3.0, -0.2, -3.0, -0.15])
59+
Q_ub = np.array([3.0, -0.2, 3.0, -0.15])
60+
theta_lb = np.array([0.0, -np.pi / 2, -np.pi / 2, -np.pi / 2])
61+
theta_ub = np.array([0.0, np.pi / 2, np.pi / 2, np.pi / 2])
62+
63+
with gp.Env() as env, gp.Model("OptimalPowerFlow", env=env) as model:
64+
P = model.addMVar(N, name="P", lb=P_lb, ub=P_ub) # Real power for buses
65+
Q = model.addMVar(N, name="Q", lb=Q_lb, ub=Q_ub) # Reactive for buses
66+
v = model.addMVar(N, name="v", lb=v_lb, ub=v_ub) # Voltage magnitude at buses
67+
68+
# Voltage angle at buses. The MVar is reshaped to a column vector to
69+
# simplify the outer subtraction used in power balance constraints.
70+
theta = model.addMVar(N, name="theta", lb=theta_lb, ub=theta_ub).reshape((N, 1))
71+
72+
# Minimize Reactive Power at buses 1 and 3
73+
model.setObjective(Q[[0, 2]].sum(), GRB.MINIMIZE)
74+
75+
# Real power balance
76+
constr_P = model.addGenConstrNL(
77+
P,
78+
v * ((G * nlfunc.cos(theta - theta.T) + B * nlfunc.sin(theta - theta.T)) @ v),
79+
name="constr_P",
80+
)
81+
82+
# Reactive power balance
83+
constr_Q = model.addGenConstrNL(
84+
Q,
85+
v * ((G * nlfunc.sin(theta - theta.T) - B * nlfunc.cos(theta - theta.T)) @ v),
86+
name="constr_Q",
87+
)
88+
89+
model.optimize()
90+
91+
# Print output table if pandas is installed
92+
try:
93+
import pandas as pd
94+
95+
df = pd.DataFrame(
96+
{
97+
"Bus": range(1, N + 1),
98+
"Voltage": v.X,
99+
"Angle": theta.reshape(N).X * 180 / np.pi,
100+
"Real Power": P.X,
101+
"Reactive Power": Q.X,
102+
}
103+
)
104+
105+
print(df)
106+
except ImportError:
107+
pass

‎cvxpy/sandbox/pypower-acopf.py‎

Lines changed: 75 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,75 @@
1+
import cvxpy as cp
2+
import numpy as np
3+
4+
N = 4
5+
6+
# Conductance/susceptance components
7+
G = np.array(
8+
[
9+
[1.7647, -0.5882, 0.0, -1.1765],
10+
[-0.5882, 1.5611, -0.3846, -0.5882],
11+
[0.0, -0.3846, 1.5611, -1.1765],
12+
[-1.1765, -0.5882, -1.1765, 2.9412],
13+
]
14+
)
15+
16+
B = np.array(
17+
[
18+
[-7.0588, 2.3529, 0.0, 4.7059],
19+
[2.3529, -6.629, 1.9231, 2.3529],
20+
[0.0, 1.9231, -6.629, 4.7059],
21+
[4.7059, 2.3529, 4.7059, -11.7647],
22+
]
23+
)
24+
25+
# Define bounds
26+
v_lb = np.array([1.0, 0.0, 1.0, 0.0])
27+
v_ub = np.array([1.0, 1.5, 1.0, 1.5])
28+
29+
P_lb = np.array([-3.0, -0.3, 0.3, -0.2])
30+
P_ub = np.array([3.0, -0.3, 0.3, -0.2])
31+
32+
Q_lb = np.array([-3.0, -0.2, -3.0, -0.15])
33+
Q_ub = np.array([3.0, -0.2, 3.0, -0.15])
34+
35+
theta_lb = np.array([0.0, -np.pi / 2, -np.pi / 2, -np.pi / 2])
36+
theta_ub = np.array([0.0, np.pi / 2, np.pi / 2, np.pi / 2])
37+
38+
# Create variables with bounds included
39+
P = cp.Variable(N, name="P", bounds=[P_lb, P_ub])
40+
Q = cp.Variable(N, name="Q", bounds=[Q_lb, Q_ub])
41+
v = cp.Variable(N, name="v", bounds=[v_lb, v_ub])
42+
theta = cp.Variable(N, name="theta", bounds=[theta_lb, theta_ub])
43+
44+
# Reshape theta to column vector for broadcasting
45+
theta_col = cp.reshape(theta, (N, 1))
46+
47+
# Create constraints list (only power balance constraints now)
48+
constraints = []
49+
50+
# Real power balance
51+
P_balance = cp.multiply(
52+
v,
53+
(
54+
G @ cp.cos(theta_col - theta_col.T)
55+
+ B @ cp.sin(theta_col - theta_col.T)
56+
) @ v
57+
)
58+
constraints.append(P == P_balance)
59+
60+
# Reactive power balance
61+
Q_balance = cp.multiply(
62+
v,
63+
(
64+
G @ cp.sin(theta_col - theta_col.T)
65+
- B @ cp.cos(theta_col - theta_col.T)
66+
) @ v
67+
)
68+
constraints.append(Q == Q_balance)
69+
70+
# Objective: minimize reactive power at buses 1 and 3 (indices 0 and 2)
71+
objective = cp.Minimize(Q[0] + Q[2])
72+
73+
# Create and solve the problem
74+
problem = cp.Problem(objective, constraints)
75+
problem.solve(solver=cp.IPOPT, verbose=True, nlp=True)

‎cvxpy/settings.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -194,7 +194,7 @@
194194
# Numerical tolerances
195195
EIGVAL_TOL = 1e-10
196196
PSD_NSD_PROJECTION_TOL = 1e-8
197-
GENERAL_PROJECTION_TOL = 1e-10
197+
GENERAL_PROJECTION_TOL = 1e-8
198198
SPARSE_PROJECTION_TOL = 1e-10
199199
ATOM_EVAL_TOL = 1e-4
200200
CHOL_SYM_TOL = 1e-14

0 commit comments

Comments
 (0)