Python cvxpy.OPTIMAL Examples

def polish(orig_prob, polish_depth=5, polish_func=None, *args, **kwargs):
    # Fix noncvx variables and solve.
    for var in get_noncvx_vars(orig_prob):
        var.value = var.z.value
    old_val = None
    for t in range(polish_depth):
        fix_constr = []
        for var in get_noncvx_vars(orig_prob):
            fix_constr += var.restrict(var.value)
        polish_prob = cvx.Problem(orig_prob.objective, orig_prob.constraints + fix_constr)
        polish_prob.solve(*args, **kwargs)
        if polish_prob.status in [cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE] and \
        (old_val is None or (old_val - polish_prob.value)/(old_val + 1) > 1e-3):
            old_val = polish_prob.value
        else:
            break

    return polish_prob.value, polish_prob.status

# Add admm method to cvx Problem. 

def solve_spectral(prob, *args, **kwargs):
    """Solve the spectral relaxation with lambda = 1.
    """

    # TODO: do this efficiently without SDP lifting

    # lifted variables and semidefinite constraint
    X = cvx.Semidef(prob.n + 1)

    W = prob.f0.homogeneous_form()
    rel_obj = cvx.Minimize(cvx.sum_entries(cvx.mul_elemwise(W, X)))

    W1 = sum([f.homogeneous_form() for f in prob.fs if f.relop == '<='])
    W2 = sum([f.homogeneous_form() for f in prob.fs if f.relop == '=='])

    rel_prob = cvx.Problem(
        rel_obj,
        [
            cvx.sum_entries(cvx.mul_elemwise(W1, X)) <= 0,
            cvx.sum_entries(cvx.mul_elemwise(W2, X)) == 0,
            X[-1, -1] == 1
        ]
    )
    rel_prob.solve(*args, **kwargs)

    if rel_prob.status not in [cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE]:
        raise Exception("Relaxation problem status: %s" % rel_prob.status)

    (w, v) = LA.eig(X.value)
    return np.sqrt(np.max(w))*np.asarray(v[:-1, np.argmax(w)]).flatten(), rel_prob.value 

def solve_sdr(prob, *args, **kwargs):
    """Solve the SDP relaxation.
    """

    # lifted variables and semidefinite constraint
    X = cvx.Semidef(prob.n + 1)

    W = prob.f0.homogeneous_form()
    rel_obj = cvx.Minimize(cvx.sum_entries(cvx.mul_elemwise(W, X)))
    rel_constr = [X[-1, -1] == 1]

    for f in prob.fs:
        W = f.homogeneous_form()
        lhs = cvx.sum_entries(cvx.mul_elemwise(W, X))
        if f.relop == '==':
            rel_constr.append(lhs == 0)
        else:
            rel_constr.append(lhs <= 0)

    rel_prob = cvx.Problem(rel_obj, rel_constr)
    rel_prob.solve(*args, **kwargs)

    if rel_prob.status not in [cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE]:
        raise Exception("Relaxation problem status: %s" % rel_prob.status)

    return X.value, rel_prob.value


# phase 1: optimize infeasibility 

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 _run_cvx_optimization(self, next_states, rewards, **solver_options):
        """Tensorflow wrapper around a cvxpy value function optimization.

        Parameters
        ----------
        next_states : ndarray
        rewards : ndarray

        Returns
        -------
        values : ndarray
            The optimal values at the states.
        """
        # Define random variables; convert index from np.int64 to regular
        # python int to avoid strange cvxpy error; see:
        # https://github.com/cvxgrp/cvxpy/issues/380
        values = cvxpy.Variable(rewards.shape)

        value_matrix = self.value_function.tri.parameter_derivative(
            next_states)
        # Make cvxpy work with sparse matrices
        value_matrix = cvxpy.Constant(value_matrix)

        objective = cvxpy.Maximize(cvxpy.sum(values))
        constraints = [values <= rewards + self.gamma * value_matrix * values]
        prob = cvxpy.Problem(objective, constraints)

        # Solve optimization problem
        prob.solve(**solver_options)

        # Some error checking
        if not prob.status == cvxpy.OPTIMAL:
            raise OptimizationError('Optimization problem is {}'
                                    .format(prob.status))

        return np.array(values.value) 

def summary(self):
        """Summary of results. Only works for single period optimization.

        :rtype: str
        """

        retval = "Status: " + self.status if self.status else "none"
        if not self.status in {cvx.OPTIMAL, cvx.OPTIMAL_INACCURATE}:
            return retval

        retval += "\n"
        retval += "%-20s %10s\n" % ("Terminal", "Power")
        retval += "%-20s %10s\n" % ("--------", "-----")
        averages = False
        for device_terminal, value in self.power.items():
            label = "%s[%d]" % (device_terminal[0].name, device_terminal[1])
            if isinstance(value, np.ndarray):
                value = np.mean(value)
                averages = True
            retval += "%-20s %10.2f\n" % (label, value)

        retval += "\n"
        retval += "%-20s %10s\n" % ("Net", "Price")
        retval += "%-20s %10s\n" % ("---", "-----")
        for net, value in self.price.items():
            if isinstance(value, np.ndarray):
                value = np.mean(value)
            retval += "%-20s %10.4f\n" % (net.name, value)

        retval += "\n"
        retval += "%-20s %10s\n" % ("Device", "Payment")
        retval += "%-20s %10s\n" % ("------", "-------")
        device_payments = {d[0][0]: 0 for d in self.payments.items()}
        for device_terminal, value in self.payments.items():
            if isinstance(value, np.ndarray):
                value = np.sum(value)
            device_payments[device_terminal[0]] += value
        for d in device_payments.keys():
            retval += "%-20s %10.2f\n" % (d.name, device_payments[d])

        if averages:
            retval += "\nPower and price are averages over the time horizon. Payment is total.\n"

        return retval 

def run_mpc(device, time_steps, predict, execute, **kwargs):
    """Execute model predictive control.

    This method executes the model predictive control loop, roughly:

    .. code:: python

      for t in time_steps:
        predict(t)
        device.problem.solve()
        execute(t)
    ..

    It is the responsibility of the provided `predict` and `execute` functions
    to update the device models with the desired predictions and execute the
    actions as appropriate.

    :param device: Device (or network of devices) to optimize
    :param time_steps: Time steps to optimize over
    :param predict: Prediction step
    :param execute: Execution step
    :type device: :class:`Device`
    :type time_steps: sequence
    :type predict: single argument function
    :type execute: single argument function
    :returns: Model predictive control results
    :rtype: :class:`Results`
    :raise: :class:`OptimizationError`

    """
    total_cost = 0.
    results = Results()
    #T_MPC = device
    for t in tqdm.trange(time_steps):
        predict(t)

        # device.init_problem(time_horizon=1)
        device.problem.solve(**kwargs)
        if device.problem.status != cvx.OPTIMAL:
            # temporary
            raise OptimizationError(
                "failed at iteration %d, %s" % (t, device.problem.status))
        stage_cost = sum([device.cost[0, 0]
                          for device in device.devices]).value
        #print('at time %s, adding cost %f' % (t, stage_cost))
        total_cost += stage_cost
        execute(t)
        _update_mpc_results(t, time_steps, device.results, results)
    return total_cost, results