add example from jing (multi-layer perceptron)

cvxpy/reductions/expr2smooth/canonicalizers/maximum_canon.py

@@ -14,10 +14,11 @@
 limitations under the License.
 """
 
 import operator
 from functools import reduce
 
+from cvxpy.atoms.affine.binary_operators import multiply
 from cvxpy.atoms.affine.wraps import nonneg_wrap, nonpos_wrap
 from cvxpy.expressions.variable import Variable
 
 
@@ -25,12 +26,8 @@
     shape = expr.shape
     t = Variable(shape)
 
-    if expr.is_nonneg():
-        t = nonneg_wrap(t)
-    if expr.is_nonpos():
-        t = nonpos_wrap(t)
-
     constraints = [t >= elem for elem in args]
     terms = [t - elem for elem in args]
-    constraints.append(reduce(operator.mul, terms) == 0)
+    expr = multiply(*terms)
+    constraints.append(expr == 0)
     return t, constraints

cvxpy/sandbox/acopf.py

@@ -0,0 +1,77 @@
+import cvxpy as cp
+import numpy as np
+
+N = 4
+
+# Conductance/susceptance components
+G = np.array(
+    [
+        [1.7647, -0.5882, 0.0, -1.1765],
+        [-0.5882, 1.5611, -0.3846, -0.5882],
+        [0.0, -0.3846, 1.5611, -1.1765],
+        [-1.1765, -0.5882, -1.1765, 2.9412],
+    ]
+)
+
+B = np.array(
+    [
+        [-7.0588, 2.3529, 0.0, 4.7059],
+        [2.3529, -6.629, 1.9231, 2.3529],
+        [0.0, 1.9231, -6.629, 4.7059],
+        [4.7059, 2.3529, 4.7059, -11.7647],
+    ]
+)
+
+# Assign bounds where fixings are needed
+v_lb = np.array([1.0, 0.0, 1.0, 0.0])
+v_ub = np.array([1.0, 1.5, 1.0, 1.5])
+
+P_lb = np.array([-3.0, -0.3, 0.3, -0.2])
+P_ub = np.array([3.0, -0.3, 0.3, -0.2])
+
+Q_lb = np.array([-3.0, -0.2, -3.0, -0.15])
+Q_ub = np.array([3.0, -0.2, 3.0, -0.15])
+
+theta_lb = np.array([0.0, -np.pi / 2, -np.pi / 2, -np.pi / 2])
+theta_ub = np.array([0.0, np.pi / 2, np.pi / 2, np.pi / 2])
+
+# Create variables with bounds
+P = cp.Variable(N, name="P")
+P.value = np.random.uniform(P_lb, P_ub)  # Real power for buses
+Q = cp.Variable(N, name="Q")
+Q.value = np.random.uniform(Q_lb, Q_ub)  # Reactive power for buses
+v = cp.Variable(N, name="v")
+v.value = np.random.uniform(v_lb, v_ub)  # Voltage magnitude at buses
+theta = cp.Variable(N, name="theta")
+theta.value = np.random.uniform(theta_lb, theta_ub)  # Voltage angle at buses
+
+# Reshape theta to column vector for broadcasting
+theta_col = cp.reshape(theta, (N, 1))
+
+# Create constraints list
+constraints = []
+
+# Add bound constraints
+constraints += [
+    P >= P_lb,
+    P <= P_ub,
+    Q >= Q_lb,
+    Q <= Q_ub,
+    v >= v_lb,
+    v <= v_ub,
+    theta >= theta_lb,
+    theta <= theta_ub
+]
+P_balance = cp.multiply(v, (G * cp.cos(theta_col - theta_col.T) + B * cp.sin(theta_col - theta_col.T)) @ v)
+constraints.append(P == P_balance)
+
+# Reactive power balance
+Q_balance = cp.multiply(v, (G * cp.sin(theta_col - theta_col.T) - B * cp.cos(theta_col - theta_col.T)) @ v)
+constraints.append(Q == Q_balance)
+
+# Objective: minimize reactive power at buses 1 and 3 (indices 0 and 2)
+objective = cp.Minimize(Q[0] + Q[2])
+
+# Create and solve the problem
+problem = cp.Problem(objective, constraints)
+problem.solve(solver=cp.IPOPT, nlp=True)

cvxpy/sandbox/standard_form.py

@@ -0,0 +1,16 @@
+import cvxpy as cp
+
+n = 5
+x = cp.Variable(n)
+# form a simple linear programming problem
+objective = cp.Minimize(cp.sum(x))
+constraints = [x >= 0, cp.sum(x) == 1]
+problem = cp.Problem(objective, constraints)
+data, chain, inverse_data = problem.get_problem_data(cp.CLARABEL)
+print(data)
+print(data['A'].toarray())
+print(data['b'])
+print("----- Chain -----")
+print(chain)
+print("----- Inverse Data -----")
+print(inverse_data)

cvxpy/sandbox/train_MLP/trainMLP.py

@@ -0,0 +1,64 @@
+import numpy as np
+import pandas as pd
+import cvxpy as cp
+import matplotlib.pyplot as plt
+
+class TrainMLP():
+    """
+    Train Multi-layer-perceptron.
+        Params:
+            [W_1, b_1, W_2, b_2]: (p + 3)-vector
+
+        Ex:
+            p = 8
+    """
+    def __init__(self, F_mat, z_vec, r_vec):
+        self.F_mat = F_mat  # T by p matrix
+        self.z_vec = z_vec  # T by 1 vector
+        self.r_vec = r_vec  # T by 1 vector
+        self.T, self.p = F_mat.shape
+
+    def loss_func(self):
+        W_1 = cp.Variable((1, self.p))  # 1 by p vector
+        b_1 = cp.Variable(1)     
+        W_2 = cp.Variable(1)    
+        b_2 = cp.Variable(1)     
+
+        MLP_vec = W_2 * cp.maximum(0, self.F_mat @ W_1.T + b_1) + b_2  # T by 1 vector
+        MLP_vec = cp.reshape(MLP_vec, (self.T,))
+        objective = cp.Minimize(cp.sum(cp.multiply(self.z_vec, cp.square(self.r_vec - MLP_vec))))
+        # No constraint in my problem. Consider it here for optimization purpose.
+        constraints = [MLP_vec >=0 , MLP_vec <= 1]
+        problem = cp.Problem(objective, constraints)
+        problem.solve(solver=cp.IPOPT, nlp=True)
+
+        return W_1.value, b_1.value, W_2.value, b_2.value
+
+    def predict(self, F_new, W_1, b_1, W_2, b_2):
+        hidden = np.maximum(0, F_new @ W_1.T + b_1)
+        predictions = W_2 * hidden + b_2
+        return predictions
+
+
+F_mat = pd.read_pickle("/Users/willizz/Documents/cvxpy-ipopt/cvxpy/sandbox/train_MLP/F_mat.pkl")
+z_vec = pd.read_pickle("/Users/willizz/Documents/cvxpy-ipopt/cvxpy/sandbox/train_MLP/z_vec.pkl")
+r_vec_1 = pd.read_pickle("/Users/willizz/Documents/cvxpy-ipopt/cvxpy/sandbox/train_MLP/r_vec_1.pkl")
+r_vec_3 = pd.read_pickle("/Users/willizz/Documents/cvxpy-ipopt/cvxpy/sandbox/train_MLP/r_vec_3.pkl")
+r_vec_5 = pd.read_pickle("/Users/willizz/Documents/cvxpy-ipopt/cvxpy/sandbox/train_MLP/r_vec_5.pkl")
+r_vec_10 = pd.read_pickle("/Users/willizz/Documents/cvxpy-ipopt/cvxpy/sandbox/train_MLP/r_vec_10.pkl")
+r_vec_50 = pd.read_pickle("/Users/willizz/Documents/cvxpy-ipopt/cvxpy/sandbox/train_MLP/r_vec_50.pkl")
+
+F_mat = F_mat.to_numpy()
+z_vec = z_vec.to_numpy().flatten()
+r_vec_1 = r_vec_1.to_numpy().flatten()
+r_vec_3 = r_vec_3.to_numpy().flatten()
+r_vec_5 = r_vec_5.to_numpy().flatten()
+r_vec_10 = r_vec_10.to_numpy().flatten()
+r_vec_50 = r_vec_50.to_numpy().flatten()
+
+MLP = TrainMLP(F_mat, z_vec, r_vec_1)
+W_1, b_1, W_2, b_2 = MLP.loss_func()
+pred = MLP.predict(F_mat, W_1, b_1, W_2, b_2)
+#print(pred[:20])
+print("W_1, b_1, W_2, b_2:")
+print(W_1, b_1, W_2, b_2)