#!/usr/bin/env python
# coding: utf-8

# # Second-order cone program

# A second-order cone program (SOCP) is an optimization problem of the form
# $$  
#     \begin{array}{ll}
#     \mbox{minimize}   & f^Tx\\
#     \mbox{subject to} & \|A_ix + b_i\|_2 \leq c_i^Tx + d_i, \quad i=1,\ldots,m \\
#                       & Fx = g,
#     \end{array}
# $$
# where $x \in \mathcal{R}^{n}$ is the optimization variable and $f \in \mathcal{R}^n$, $A_i \in \mathcal{R}^{n_i \times n}$, $b_i \in \mathcal{R}^{n_i}$, $c_i \in \mathcal{R}^n$, $d_i \in \mathcal{R}$, $F \in \mathcal{R}^{p \times n}$, and $g \in \mathcal{R}^p$ are problem data.
# 
# An example of an SOCP is the robust linear program
# $$  
#     \begin{array}{ll}
#     \mbox{minimize}   & c^Tx\\
#     \mbox{subject to} & (a_i + u_i)^Tx \leq b_i \textrm{ for all } \|u_i\|_2 \leq 1, \quad i=1,\ldots,m,
#     \end{array}
# $$
# where the problem data $a_i$ are known within an $\ell_2$-norm ball of radius one. The robust linear program can be rewritten as the SOCP
# $$  
#     \begin{array}{ll}
#     \mbox{minimize}   & c^Tx\\
#     \mbox{subject to} & a_i^Tx + \|x\|_2 \leq b_i, \quad i=1,\ldots,m,
#     \end{array}
# $$
# 
# When we solve a SOCP, in addition to a solution $x^\star$, we obtain a dual solution $\lambda_i^\star$ corresponding to each second-order cone constraint. A non-zero $\lambda_i^\star$ indicates that the constraint $ \|A_ix + b_i\|_2 \leq c_i^Tx + d_i$ holds with equality for $x^\star$ and suggests that changing $d_i$ would change the optimal value.

# ## Example

# In the following code, we solve an SOCP with CVXPY.

# In[31]:


# Import packages.
import cvxpy as cp
import numpy as np

# Generate a random feasible SOCP.
m = 3
n = 10
p = 5
n_i = 5
np.random.seed(2)
f = np.random.randn(n)
A = []
b = []
c = []
d = []
x0 = np.random.randn(n)
for i in range(m):
    A.append(np.random.randn(n_i, n))
    b.append(np.random.randn(n_i))
    c.append(np.random.randn(n))
    d.append(np.linalg.norm(A[i]@x0 + b, 2) - c[i].T@x0)
F = np.random.randn(p, n)
g = F@x0

# Define and solve the CVXPY problem.
x = cp.Variable(n)
# We use cp.SOC(t, x) to create the SOC constraint ||x||_2 <= t.
soc_constraints = [
    cp.SOC(c[i].T@x + d[i], A[i]@x + b[i]) for i in range(m)
]
prob = cp.Problem(cp.Minimize(f.T@x),
                  soc_constraints + [F@x == g])
prob.solve()

# Print result.
print("The optimal value is", prob.value)
print("A solution x is")
print(x.value)
for i in range(m):
    print("SOC constraint %i dual variable solution" % i)
    print(soc_constraints[i].dual_value)