Python cvxpy.quad_form() Examples

def _generate_cvxpy_problem(self):
        '''
        Generate QP problem
        '''

        n = self.n
        m = self.m
        x = cvxpy.Variable(n)
        t = cvxpy.Variable(m)

        objective = cvxpy.Minimize(.5 * cvxpy.quad_form(x, spa.eye(n))
                                   + .5 * self.gamma * np.ones(m) * t)
        constraints = [t >= spa.diags(self.b_svm).dot(self.A_svm) * x + 1,
                       t >= 0]

        problem = cvxpy.Problem(objective, constraints)

        return problem, (x, t) 

def _generate_cvxpy_problem(self):
        '''
        Generate QP problem
        '''

        x = cvxpy.Variable(self.n)
        y = cvxpy.Variable(self.k)

        # Create parameters m
        mu = cvxpy.Parameter(self.n)
        mu.value = self.mu

        objective = cvxpy.Minimize(cvxpy.quad_form(x, self.D) +
                                   cvxpy.quad_form(y, spa.eye(self.k)) +
                                   - 1 / self.gamma * (mu.T * x))
        constraints = [np.ones(self.n) * x == 1,
                       self.F.T * x == y,
                       0 <= x, x <= 1]
        problem = cvxpy.Problem(objective, constraints)

        return problem, mu 

def eval_cvx(self, x):
        return cvx.quad_form(x, self.P) + self.q.T*x + self.r 

def forward_single_np(Q, p, G, h, A, b):
    nz, neq, nineq = p.shape[0], A.shape[0] if A is not None else 0, G.shape[0]

    z_ = cp.Variable(nz)

    obj = cp.Minimize(0.5 * cp.quad_form(z_, Q) + p.T * z_)
    eqCon = A * z_ == b if neq > 0 else None
    if nineq > 0:
        slacks = cp.Variable(nineq)
        ineqCon = G * z_ + slacks == h
        slacksCon = slacks >= 0
    else:
        ineqCon = slacks = slacksCon = None
    cons = [x for x in [eqCon, ineqCon, slacksCon] if x is not None]
    prob = cp.Problem(obj, cons)
    prob.solve()  # solver=cp.SCS, max_iters=5000, verbose=False)
    # prob.solve(solver=cp.SCS, max_iters=10000, verbose=True)
    assert('optimal' in prob.status)
    zhat = np.array(z_.value).ravel()
    nu = np.array(eqCon.dual_value).ravel() if eqCon is not None else None
    if ineqCon is not None:
        lam = np.array(ineqCon.dual_value).ravel()
        slacks = np.array(slacks.value).ravel()
    else:
        lam = slacks = None

    return prob.value, zhat, nu, lam, slacks 

def mpc_control(x0):

    x = cvxpy.Variable((nx, T + 1))
    u = cvxpy.Variable((nu, T))

    A, B = get_model_matrix()

    cost = 0.0
    constr = []
    for t in range(T):
        cost += cvxpy.quad_form(x[:, t + 1], Q)
        cost += cvxpy.quad_form(u[:, t], R)
        constr += [x[:, t + 1] == A * x[:, t] + B * u[:, t]]
    # print(x0)
    constr += [x[:, 0] == x0[:, 0]]
    prob = cvxpy.Problem(cvxpy.Minimize(cost), constr)

    start = time.time()
    prob.solve(verbose=False)
    elapsed_time = time.time() - start
    print("calc time:{0} [sec]".format(elapsed_time))

    if prob.status == cvxpy.OPTIMAL:
        ox = get_nparray_from_matrix(x.value[0, :])
        dx = get_nparray_from_matrix(x.value[1, :])
        theta = get_nparray_from_matrix(x.value[2, :])
        dtheta = get_nparray_from_matrix(x.value[3, :])

        ou = get_nparray_from_matrix(u.value[0, :])

    return ox, dx, theta, dtheta, ou 

def use_modeling_tool(A, B, N, Q, R, P, x0, umax=None, umin=None, xmin=None, xmax=None):
    """
    solve MPC with modeling tool for test
    """
    (nx, nu) = B.shape

    # mpc calculation
    x = cvxpy.Variable((nx, N + 1))
    u = cvxpy.Variable((nu, N))

    costlist = 0.0
    constrlist = []

    for t in range(N):
        costlist += 0.5 * cvxpy.quad_form(x[:, t], Q)
        costlist += 0.5 * cvxpy.quad_form(u[:, t], R)

        constrlist += [x[:, t + 1] == A * x[:, t] + B * u[:, t]]

        if xmin is not None:
            constrlist += [x[:, t] >= xmin[:, 0]]
        if xmax is not None:
            constrlist += [x[:, t] <= xmax[:, 0]]

    costlist += 0.5 * cvxpy.quad_form(x[:, N], P)  # terminal cost
    if xmin is not None:
        constrlist += [x[:, N] >= xmin[:, 0]]
    if xmax is not None:
        constrlist += [x[:, N] <= xmax[:, 0]]

    constrlist += [x[:, 0] == x0[:, 0]]  # inital state constraints
    if umax is not None:
        constrlist += [u <= umax]  # input constraints
    if umin is not None:
        constrlist += [u >= umin]  # input constraints

    prob = cvxpy.Problem(cvxpy.Minimize(costlist), constrlist)

    prob.solve(verbose=True)

    return x.value, u.value 

def full_qp():
    # print(f'--- {sys._getframe().f_code.co_name} ---')
    print('full qp')
    npr.seed(0)
    nx, ncon_eq, ncon_ineq = 5, 2, 3

    Q = cp.Parameter((nx, nx))
    p = cp.Parameter((nx, 1))
    A = cp.Parameter((ncon_eq, nx))
    b = cp.Parameter(ncon_eq)
    G = cp.Parameter((ncon_ineq, nx))
    h = cp.Parameter(ncon_ineq)
    x = cp.Variable(nx)
    # obj = cp.Minimize(0.5*cp.quad_form(x, Q) + p.T * x)
    obj = cp.Minimize(0.5 * cp.sum_squares(Q@x) + p.T * x)
    cons = [A * x == b, G * x <= h]
    prob = cp.Problem(obj, cons)

    x0 = npr.randn(nx)
    s0 = npr.randn(ncon_ineq)

    G.value = npr.randn(ncon_ineq, nx)
    h.value = G.value.dot(x0) + s0

    A.value = npr.randn(ncon_eq, nx)
    b.value = A.value.dot(x0)

    L = npr.randn(nx, nx)
    Q.value = L.T  # L.dot(L.T)
    p.value = npr.randn(nx, 1)

    prob.solve(solver=cp.SCS)
    print(x.value) 

def _generate_cvxpy_problem(self):
        '''
        Generate QP problem
        '''
        x_var = cvxpy.Variable(self.n)
        objective = .5 * cvxpy.quad_form(x_var, self.P) + self.q * x_var
        constraints = [self.A * x_var == self.u]
        problem = cvxpy.Problem(cvxpy.Minimize(objective), constraints)

        return problem 

def _generate_cvxpy_problem(self):
        '''
        Generate QP problem
        '''
        x_var = cvxpy.Variable(self.n)
        objective = .5 * cvxpy.quad_form(x_var, self.P) + self.q * x_var + \
            self.r
        constraints = [self.A * x_var <= self.u, self.A * x_var >= self.l]
        problem = cvxpy.Problem(cvxpy.Minimize(objective), constraints)

        return problem 

def _generate_cvxpy_problem(self):
        '''
        Generate QP problem
        '''
        x_var = cvxpy.Variable(self.n)
        objective = .5 * cvxpy.quad_form(x_var, self.P) + self.q * x_var + \
            self.r
        constraints = [self.A * x_var <= self.u, self.A * x_var >= self.l]
        problem = cvxpy.Problem(cvxpy.Minimize(objective), constraints)

        return problem 

def _generate_cvxpy_problem(self):
        '''
        Generate QP problem
        '''
        x_var = cvxpy.Variable(self.n)
        objective = .5 * cvxpy.quad_form(x_var, self.P) + self.q * x_var
        constraints = [self.A * x_var <= self.u, self.A * x_var >= self.l]
        problem = cvxpy.Problem(cvxpy.Minimize(objective), constraints)

        return problem 

def portfolio_opt(p):
    """ Create a portfolio optimization problem with p dimensions """
    temp = np.random.randn(p, p)
    Sigma = temp.T.dot(temp)


    Sigma_sqrt = sp.linalg.sqrtm(Sigma)
    o = np.ones((p, 1))

    # Create standard form cone problem
    Zp1 = np.zeros((p,1))
    
    # setup for cone problem
    c = np.vstack([Zp1, np.array([[1.0]])]).ravel()

    G1 = sp.linalg.block_diag(2.0*Sigma_sqrt, -1.0)
    q = np.vstack([Zp1, np.array([[1.0]])])
    G2 = np.hstack([o.T, np.array([[0.0]])])
    G3 = np.hstack([-o.T, np.array([[0.0]])])

    h = np.vstack([Zp1, np.array([[1.0]])])
    z = 1.0

    A = np.vstack([G2, G3, -q.T, -G1 ])
    b = np.vstack([1.0, -1.0, np.array([[z]]), h]).ravel()

    betahat = cp.Variable(p)
    return (betahat, cp.Problem(cp.Minimize(cp.quad_form(betahat, Sigma)),
                                [o.T * betahat == 1]), {
        'A' : A, 
        'b' : b,
        'c' : c, 
        'dims' : {
            'l' : 2,
            'q' : [p+2]
        },
        'beta_from_x' : lambda x: x[:p]
    }) 

def optimize_focality_cvxpy(l, Q, target_mean, max_total_current=None,
                    max_el_current=None, Qin=None, max_angle=None, return_duals=False,
                    none_on_infeasibility=False):
    import cvxpy

    if max_total_current is None and max_el_current is None:
        raise ValueError('Please define a maximal total current or maximal electrode ' +
                         'current')
    n = l.shape[0]
    C = np.vstack([np.eye(n), -np.eye(n)])
    A = np.ones(n)

    if max_el_current is not None:
        d = max_el_current * np.ones(C.shape[0])
        eps = 1e-3 * max_el_current
    else:
        d = 1e10 * np.ones(C.shape[0])

    if max_total_current is not None:
        max_l1 = 2 * max_total_current
    else:
        max_l1 = None

    C = np.vstack([C, l])
    d = np.hstack([d, target_mean])
    x = cvxpy.Variable(n)
    p = cvxpy.Problem(cvxpy.Maximize(l * x))
    p.constraints = [C * x <= d,
                     A * x == 0]
    if max_total_current is not None:
        p.constraints.append(cvxpy.norm(x, 1) <= max_l1)
    p.solve(solver=cvxpy.SCS)
    v = np.squeeze(np.array(x.value)).T

    field_component = l.dot(v)
    if field_component * (1 + 1e-4) < target_mean:
        return v
    else:
        target_field = target_mean
        # Solve the QP
        eq_constraints = np.vstack([A, l])
        b = np.array([0, target_field])
        C = C[:-1]
        d = d[:-1]
        p = cvxpy.Problem(cvxpy.Minimize(cvxpy.quad_form(x, Q)))
        p.constraints = [C * x <= d,
                         eq_constraints * x == b]
        if max_total_current is not None:
            p.constraints.append(cvxpy.norm(x, 1) <= max_l1)
        p.solve(solver=cvxpy.SCS)
        v = np.squeeze(np.array(x.value)).T

        return v 

def cvxpy_solve_qp(P, q, G=None, h=None, A=None, b=None, initvals=None,
                   solver=None, verbose=False):
    """
    Solve a Quadratic Program defined as:

    .. math::

        \\begin{split}\\begin{array}{ll}
        \\mbox{minimize} &
            \\frac{1}{2} x^T P x + q^T x \\\\
        \\mbox{subject to}
            & G x \\leq h                \\\\
            & A x = h
        \\end{array}\\end{split}

    calling a given solver using the `CVXPY <http://www.cvxpy.org/>`_ modelling
    language.

    Parameters
    ----------
    P : array, shape=(n, n)
        Primal quadratic cost matrix.
    q : array, shape=(n,)
        Primal quadratic cost vector.
    G : array, shape=(m, n)
        Linear inequality constraint matrix.
    h : array, shape=(m,)
        Linear inequality constraint vector.
    A : array, shape=(meq, n), optional
        Linear equality constraint matrix.
    b : array, shape=(meq,), optional
        Linear equality constraint vector.
    initvals : array, shape=(n,), optional
        Warm-start guess vector (not used).
    solver : string, optional
        Solver name in ``cvxpy.installed_solvers()``.
    verbose : bool, optional
        Set to `True` to print out extra information.

    Returns
    -------
    x : array, shape=(n,)
        Solution to the QP, if found, otherwise ``None``.
    """
    if initvals is not None:
        print("CVXPY: note that warm-start values are ignored by wrapper")
    n = q.shape[0]
    x = Variable(n)
    P = Constant(P)  # see http://www.cvxpy.org/en/latest/faq/
    objective = Minimize(0.5 * quad_form(x, P) + q * x)
    constraints = []
    if G is not None:
        constraints.append(G * x <= h)
    if A is not None:
        constraints.append(A * x == b)
    prob = Problem(objective, constraints)
    prob.solve(solver=solver, verbose=verbose)
    x_opt = array(x.value).reshape((n,))
    return x_opt 

def solve_quadratic_program(self, P, q):
        """
        Solve the quadratic program optimization problem.

        This function solved the quadratic program to minimize the objective function
        f(x) = 1/2(x*P*x)+q*x
        subject to the additional constraints
        Gx <= h

        Where P, q are given and G,h are computed to ensure that x represents
        a probability vector and subject to honesty constraints if required
        Args:
            P (matrix): A matrix representing the P component of the objective function
            q (list): A vector representing the q component of the objective function

        Returns:
            list: The solution of the quadratic program (represents probabilities)

        Additional information:
            This method is the only place in the code where we rely on the cvxpy library
            should we consider another library, only this method needs to change.
        """
        try:
            import cvxpy
        except ImportError:
            logger.error("cvxpy module needs to be installed to use this feature.")

        P = numpy.array(P).astype(float)
        q = numpy.array(q).astype(float).T
        n = len(q)
        # G and h constrain:
        #   1) sum of probs is less then 1
        #   2) All probs bigger than 0
        #   3) Honesty (measured using fidelity, if given)
        G_data = [[1] * n] + [([-1 if i == k else 0 for i in range(n)])
                              for k in range(n)]
        h_data = [1] + [0] * n
        if self.fidelity_data is not None:
            G_data.append(self.fidelity_data['coefficients'])
            h_data.append(self.fidelity_data['goal'])
        G = numpy.array(G_data).astype(float)
        h = numpy.array(h_data).astype(float)
        x = cvxpy.Variable(n)
        prob = cvxpy.Problem(
            cvxpy.Minimize((1 / 2) * cvxpy.quad_form(x, P) + q.T @ x),
            [G @ x <= h])
        prob.solve()
        return x.value 

def use_modeling_tool(A, B, N, Q, R, P, x0, umax=None, umin=None, xmin=None, xmax=None):
    """
    solve MPC with modeling tool for test
    """
    (nx, nu) = B.shape

    # mpc calculation
    x = cvxpy.Variable((nx, N + 1))
    u = cvxpy.Variable((nu, N))

    costlist = 0.0
    constrlist = []

    for t in range(N):
        costlist += 0.5 * cvxpy.quad_form(x[:, t], Q)
        costlist += 0.5 * cvxpy.quad_form(u[:, t], R)

        constrlist += [x[:, t + 1] == A * x[:, t] + B * u[:, t]]

        if xmin is not None:
            constrlist += [x[:, t] >= xmin[:, 0]]
        if xmax is not None:
            constrlist += [x[:, t] <= xmax[:, 0]]

    costlist += 0.5 * cvxpy.quad_form(x[:, N], P)  # terminal cost
    if xmin is not None:
        constrlist += [x[:, N] >= xmin[:, 0]]
    if xmax is not None:
        constrlist += [x[:, N] <= xmax[:, 0]]

    if umax is not None:
        constrlist += [u <= umax]  # input constraints
    if umin is not None:
        constrlist += [u >= umin]  # input constraints

    constrlist += [x[:, 0] == x0[:, 0]]  # inital state constraints

    prob = cvxpy.Problem(cvxpy.Minimize(costlist), constrlist)

    prob.solve(verbose=True)

    return x.value, u.value