Python cvxpy.Problem类代码示例

def GetMinimalSpeedToReachEpsilonNeighbordhoodVector(dt, epsilon, W, dW, dF):
        Ndim = W.shape[0]
        Nsamples = W.shape[1]

        dist_w = np.zeros((Nsamples-1))
        for i in range(0,Nsamples-1):
                dist_w[i] = np.linalg.norm(W[:,i]-W[:,i+1])

        p = Variable(Nsamples-1)
        sM = Variable(Nsamples-1)

        constraints = []
        objfunc = 0.0
        for i in range(0,Nsamples-1):
                #constraints.append( norm(p*dt*dW0[0:2] + dF +np.dot(dw,np.array((1,0))) ) < epsilon )
                constraints.append( norm(p[i]*dt*dW[:,i] + dt*dt/2*dF[:,i] + np.dot(dist_w[i],np.array((1,0,0,0))) ) < epsilon )
                constraints.append( sM[i] >= p[i] )
                constraints.append( sM[i] >= 0.0)
                constraints.append( p[i] >= 0.0 )
                objfunc += norm(sM[i])

        objective = Minimize(objfunc)

        prob = Problem(objective, constraints)

        print "solve minimal speed"
        result = prob.solve(solver=SCS)
        print "done.(",prob.value,"|",np.min(sM.value),")"

        if prob.value < inf:
                return np.array(sM.value).flatten()
        else:
                return np.array(sM.value).flatten()

    def test_partial_optimize_numeric_fn(self):
        x, y = Variable(1), Variable(1)
        xval = 4

        # Solve the (simple) two-stage problem by "combining" the two stages (i.e., by solving a single linear program)
        p1 = Problem(Minimize(y), [xval + y >= 3])
        p1.solve()

        # Solve the two-stage problem via partial_optimize
        constr = [y >= -100]
        p2 = Problem(Minimize(y), [x + y >= 3] + constr)
        g = cvxpy.partial_optimize(p2, [y], [x])
        x.value = xval
        y.value = 42
        constr[0].dual_variable.value = 42
        result = g.value
        self.assertAlmostEqual(result, p1.value)
        self.assertAlmostEqual(y.value, 42)
        self.assertAlmostEqual(constr[0].dual_value, 42)

        # No variables optimized over.
        p2 = Problem(Minimize(y), [x + y >= 3])
        g = cvxpy.partial_optimize(p2, [], [x, y])
        x.value = xval
        y.value = 42
        p2.constraints[0].dual_variable.value = 42
        result = g.value
        self.assertAlmostEqual(result, y.value)
        self.assertAlmostEqual(y.value, 42)
        self.assertAlmostEqual(p2.constraints[0].dual_value, 42)

 def test_pnorm(self):
     """ Test domain for pnorm.
     """
     dom = pnorm(self.a, -0.5).domain
     prob = Problem(Minimize(self.a), dom)
     prob.solve()
     self.assertAlmostEqual(prob.value, 0)

 def test_indicator(self):
     """Test indicator transform.
     """
     cons = [self.a >= 0, self.x == 2]
     obj = cvxpy.Minimize(self.a - sum_entries(self.x))
     expr = cvxpy.indicator(cons)
     assert expr.is_convex()
     assert expr.is_positive()
     prob = Problem(Minimize(expr) + obj)
     result = prob.solve()
     self.assertAlmostEqual(-4, result)
     self.assertAlmostEqual(0, self.a.value)
     self.assertItemsAlmostEqual([2, 2], self.x.value)
     self.assertAlmostEqual(0, expr.value)
     self.a.value = -1
     self.assertAlmostEqual(np.infty, expr.value)

 def test_matrix_frac(self):
     """Test domain for matrix_frac.
     """
     dom = matrix_frac(self.x, self.A + np.eye(2)).domain
     prob = Problem(Minimize(sum_entries(diag(self.A))), dom)
     prob.solve(solver=cvxpy.SCS)
     self.assertAlmostEquals(prob.value, -2, places=3)

    def test_partial_problem(self):
        """Test domain for partial minimization/maximization problems.
        """
        for obj in [Minimize((self.a)**-1), Maximize(log(self.a))]:
            prob = Problem(obj, [self.x + self.a >= [5, 8]])
            # Optimize over nothing.
            expr = cvxpy.partial_optimize(prob, dont_opt_vars=[self.x, self.a])
            dom = expr.domain
            constr = [self.a >= -100, self.x >= 0]
            prob = Problem(Minimize(sum_entries(self.x + self.a)), dom + constr)
            prob.solve()
            self.assertAlmostEqual(prob.value, 13)
            assert self.a.value >= 0
            assert np.all((self.x + self.a - [5, 8]).value >= -1e-3)

            # Optimize over x.
            expr = cvxpy.partial_optimize(prob, opt_vars=[self.x])
            dom = expr.domain
            constr = [self.a >= -100, self.x >= 0]
            prob = Problem(Minimize(sum_entries(self.x + self.a)), dom + constr)
            prob.solve()
            self.assertAlmostEqual(prob.value, 0)
            assert self.a.value >= 0
            self.assertItemsAlmostEqual(self.x.value, [0, 0])

            # Optimize over x and a.
            expr = cvxpy.partial_optimize(prob, opt_vars=[self.x, self.a])
            dom = expr.domain
            constr = [self.a >= -100, self.x >= 0]
            prob = Problem(Minimize(sum_entries(self.x + self.a)), dom + constr)
            prob.solve()
            self.assertAlmostEqual(self.a.value, -100)
            self.assertItemsAlmostEqual(self.x.value, [0, 0])

def ForceCanBeCounteractedByAccelerationVector(dt, Fp, u1min, u1max, u2min, u2max, plot=False) :

        ### question (1) : can we produce an acceleration ddW, such that it counteracts F?

        ## dynamics projected onto identity element, it becomes obvious that in an infinitesimal neighborhood, 
        ## we can only counteract forces along the x and the theta axes due to non-holonomicity

        dt2 = dt*dt/2

        ## span dt2-hyperball in Ndim
        F = dt2*Fp
        thetamin = dt2*u2min
        thetamax = dt2*u2max
        xmin = 0.0
        xmax = dt2*u1max

        Xlow = np.dot(np.dot(Rz(-pi/2),Rz(thetamin)),np.array((1,0,0)))
        Xhigh = np.dot(np.dot(Rz(pi/2),Rz(thetamax)),np.array((1,0,0)))

        Ndim = Fp.shape[0]
        if Fp.ndim <= 1:
                Nsamples = 1
        else:
                Nsamples = Fp.shape[1]
        p = Variable(3,Nsamples)

        constraints = []
        objfunc = 0.0
        for i in range(0,Nsamples):
                constraints.append( norm(p[:,i]) <= xmax )
                constraints.append( np.matrix(Xlow[0:3])*p[:,i] <= 0 )
                constraints.append( np.matrix(Xhigh[0:3])*p[:,i] <= 0 )
                if Fp.ndim <= 1:
                        objfunc += norm(p[:,i]-F[0:3])
                else:
                        objfunc += norm(p[:,i]-F[0:3,i])
                #objfunc.append(norm(p[:,i]-F[:,i]))

        objective = Minimize(objfunc)
        prob = Problem(objective, constraints)

        result = prob.solve(solver=SCS, eps=1e-7)

        #nearest_ddq = np.array(p.value)
        nearest_ddq = np.array(p.value/dt2)

        codimension = Ndim-nearest_ddq.shape[0]

        #print Ndim, nearest_ddq.shape
        #print codimension
        zero_rows = np.zeros((codimension,Nsamples))

        if nearest_ddq.shape[0] < Ndim:
                nearest_ddq = np.vstack((nearest_ddq,zero_rows))

        if plot:
                PlotReachableSetForceDistance(dt, u1min, u1max, u2min, u2max, -F, dt2*nearest_ddq)

        return nearest_ddq

def l1_solution(A, b, lam=0.5):
    N = A.shape[0]
    x = Variable(N)
    objective = Minimize(sum_entries(square(A * x - b)) + lam * norm(x, 1))
    constraints = []
    prob = Problem(objective, constraints)

    prob.solve()
    xhat = x.value
    return xhat

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

        minimize
            (1/2) * x.T * P * x + q.T * x

        subject to
            G * x <= h
            A * x == b

    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()``.

    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)
    x_opt = array(x.value).reshape((n,))
    return x_opt

    def test_partial_optimize_numeric_fn(self):
        x,y = Variable(1), Variable(1)
        xval = 4

        # Solve the (simple) two-stage problem by "combining" the two stages (i.e., by solving a single linear program)
        p1 = Problem(Minimize(y), [xval+y>=3])
        p1.solve()

        # Solve the two-stage problem via partial_optimize
        p2 = Problem(Minimize(y), [x+y>=3])
        g = partial_optimize(p2, [y], [x])
        result = g(x).numeric([xval])
        self.assertAlmostEqual(result, p1.value)

    def test_yield_constr_cost_min(self):
        # Create problem data.
        n = 10
        c = numpy.random.randn(n)
        P, q, r = numpy.eye(n), numpy.random.randn(n), numpy.random.randn()
        mu, Sigma = numpy.zeros(n), 0.1*numpy.eye(n)
        omega = NormalRandomVariable(mu, Sigma)
        m, eta = 100, 0.95

        # Create and solve optimization problem.
        x = Variable(n)
        yield_constr = prob(quad_form(x+omega,P)
                        + (x+omega).T*q + r >= 0, m) <= 1-eta
        p = Problem(Minimize(x.T*c), [yield_constr])
        p.solve()
        self.assert_feas(p)

    def test_partial_problem(self):
        """Test grad for partial minimization/maximization problems.
        """
        for obj in [Minimize((self.a)**-1), Maximize(entr(self.a))]:
            prob = Problem(obj, [self.x + self.a >= [5, 8]])
            # Optimize over nothing.
            expr = cvxpy.partial_optimize(prob, dont_opt_vars=[self.x, self.a])
            self.a.value = None
            self.x.value = None
            grad = expr.grad
            self.assertAlmostEqual(grad[self.a], None)
            self.assertAlmostEqual(grad[self.x], None)
            # Outside domain.
            self.a.value = 1.0
            self.x.value = [5, 5]
            grad = expr.grad
            self.assertAlmostEqual(grad[self.a], None)
            self.assertAlmostEqual(grad[self.x], None)

            self.a.value = 1
            self.x.value = [10, 10]
            grad = expr.grad
            self.assertAlmostEqual(grad[self.a], obj.args[0].grad[self.a])
            self.assertItemsAlmostEqual(grad[self.x].todense(), [0, 0, 0, 0])

            # Optimize over x.
            expr = cvxpy.partial_optimize(prob, opt_vars=[self.x])
            self.a.value = 1
            grad = expr.grad
            self.assertAlmostEqual(grad[self.a], obj.args[0].grad[self.a] + 0)

            # Optimize over a.
            fix_prob = Problem(obj, [self.x + self.a >= [5, 8], self.x == 0])
            fix_prob.solve()
            dual_val = fix_prob.constraints[0].dual_variable.value
            expr = cvxpy.partial_optimize(prob, opt_vars=[self.a])
            self.x.value = [0, 0]
            grad = expr.grad
            self.assertItemsAlmostEqual(grad[self.x].todense(), dual_val)

            # Optimize over x and a.
            expr = cvxpy.partial_optimize(prob, opt_vars=[self.x, self.a])
            grad = expr.grad
            self.assertAlmostEqual(grad, {})

def mcFrobSolveLeftFactor_cvx(V, M_Omega, mask, **kwargs):
    """
    mcFrobSolveLeftFactor_cvx(V, M_Omega, mask, **kwargs)
    A solver for the left factor, U, in the problem
        min FrobNorm( P_Omega(U * V.T - M) )
    where U is an m-by-r matrix, V an n-by-r matrix.
    M_Omega is the set of observed entries in matrix form, while
    mask is a Boolean array with 1/True-valued entries corresponding 
    to those indices that were observed.

    This function is computed using the CVXPY package (and 
    thus is likely to be slower than a straight iterative 
    least squares solver).
    """
    # Options
    returnObjectiveValue = kwargs.get('returnObjectiveValue', False)
    solver = kwargs.get('solver', SCS)
    verbose = kwargs.get('verbose', False)

    if isinstance(verbose, int):
        if verbose > 1:
            verbose = True
        else:
            verbose = False

    # Parameters
    m = mask.shape[0]
    if V.shape[0] < V.shape[1]:
        # make sure V_T is "short and fat"
        V = V.T
    r = V.shape[1]

    Omega_i, Omega_j = matIndicesFromMask(mask)

    # Problem
    U = Variable(m, r)
    obj = Minimize(cvxnorm(cvxvec((U @ V.T)[Omega_i, Omega_j]) - M_Omega))
    prob = Problem(obj)
    prob.solve(solver=solver, verbose=verbose)
    if returnObjectiveValue:
        return (U.value, prob.value)
    else:
        return U.value

    def test_simple_problem(self):
        # Create problem data.
        n = numpy.random.randint(1,10)
        eta = 0.95
        num_samples = 10

        c = numpy.random.rand(n,1)

        mu = numpy.zeros(n)
        Sigma = numpy.eye(n)
        a = NormalRandomVariable(mu, Sigma)

        b = numpy.random.randn()

        # Create and solve optimization problem.
        x = Variable(n)
        p = Problem(Maximize(x.T*c), [prob(max_entries(x.T*a-b) >= 0, num_samples) <= 1-eta])
        p.solve()
        self.assert_feas(p)

    def test_value_at_risk(self):
        # Create problem data.
        n = numpy.random.randint(1,10)
        pbar = numpy.random.randn(n)
        Sigma = numpy.eye(n)
        p = NormalRandomVariable(pbar,Sigma)

        o = numpy.ones((n,1))
        beta = 0.05
        num_samples = 50

        # Create and solve optimization problem.
        x = Variable(n)
        p1 = Problem(Minimize(-x.T*pbar), [prob(-x.T*p >= 0, num_samples) <= beta, x.T*o == 1, x >= -0.1])
        p1.solve()

        # Create and solve analytic form of optimization problem (as a check).
        p2 = Problem(Minimize(-x.T*pbar),
                     [x.T*pbar >= scipy.stats.norm.ppf(1-beta) * norm2(sqrtm(Sigma) * x), x.T*o == 1, x >= -0.1])
        p2.solve()

        tol = 0.1
        if numpy.abs(p1.value - p2.value) < tol:
            self.assertAlmostEqual(1,1)
        else:
            self.assertAlmostEqual(1,0)