How to use the cvxpy.multiply function in cvxpy

s.t. x_{t+1} = A_t x_t + B_t u_t + f_t
                             x_0 = x_init
                             u_lower <= u <= u_upper
    """
    tau = cp.Variable((n_state+n_ctrl, T))
    assert (u_lower is None) == (u_upper is None)

    objs = []
    x0 = tau[:n_state,0]
    u0 = tau[n_state:,0]
    cons = [x0 == x_init]
    for t in range(T):
        xt = tau[:n_state,t]
        ut = tau[n_state:,t]
        objs.append(0.5*cp.quad_form(tau[:,t], C[t]) +
                    cp.sum(cp.multiply(c[t], tau[:,t])))
        if u_lower is not None:
            cons += [u_lower[t] <= ut, ut <= u_upper[t]]
        if t+1 < T:
            xtp1 = tau[:n_state, t+1]
            cons.append(xtp1 == F[t]*tau[:,t]+f[t])
    prob = cp.Problem(cp.Minimize(sum(objs)), cons)
    # prob.solve(solver=cp.SCS, verbose=True)
    prob.solve()
    assert 'optimal' in prob.status
    return np.array(tau.value), np.array([obj_t.value for obj_t in objs])

return 1 / (1 + np.exp(-z))

        X_np = np.random.randn(N, n)
        a_true = np.random.randn(n, 1)
        y_np = np.round(sigmoid(X_np @ a_true + np.random.randn(N, 1) * 0.5))

        X_tf = tf.Variable(X_np)
        lam_tf = tf.Variable(1.0 * tf.ones(1))

        a = cp.Variable((n, 1))
        X = cp.Parameter((N, n))
        lam = cp.Parameter(1, nonneg=True)
        y = y_np

        log_likelihood = cp.sum(
            cp.multiply(y, X @ a) -
            cp.log_sum_exp(cp.hstack([np.zeros((N, 1)), X @ a]).T, axis=0,
                           keepdims=True).T
        )
        prob = cp.Problem(
            cp.Minimize(-log_likelihood + lam * cp.sum_squares(a)))
        fit_logreg = CvxpyLayer(prob, [X, lam], [a])

        with tf.GradientTape(persistent=True) as tape:
            weights = fit_logreg(X_tf, lam_tf, solver_args={'eps': 1e-8})[0]
            summed = tf.math.reduce_sum(weights)
        grad_X_tf, grad_lam_tf = tape.gradient(summed, [X_tf, lam_tf])

        def f_train():
            prob.solve(solver=cp.SCS, eps=1e-8)
            return np.sum(a.value)

constr += [z == 0]
                second_term = 0
        else:  # it is a pd series
            no_trade = second_term.index[second_term.isnull()]
            second_term[no_trade] = 0
            constr += [z[second_term.index.get_loc(tick)] == 0
                       for tick in no_trade]

        try:
            self.expression = cvx.multiply(
                time_locator(self.half_spread, t), cvx.abs(z))
        except TypeError:
            self.expression = cvx.multiply(
                time_locator(self.half_spread, t).values, cvx.abs(z))
        try:
            self.expression += cvx.multiply(second_term,
                                            cvx.abs(z) ** self.power)
        except TypeError:
            self.expression += cvx.multiply(
                second_term.values, cvx.abs(z) ** self.power)

        return cvx.sum(self.expression), constr

for i in range(n):
            end_resource[i][tau] = np.sum(Zval[:, tau*n+i])
            lambda_right[i][1+tau] = dual_right_fut[tau*n+i][0]

    ## Generating left derivative
    small = 1.0e-6
    R_left = np.maximum(R - 0.0001 * np.ones((n, 1)), small*np.ones((n, 1)))  # R_left should be positive
    Ru_left = np.maximum(np.reshape(Ru, (tau_max*n, 1), order='F') - 0.0001 * np.ones((tau_max*n, 1)), small*np.ones((tau_max*n, 1)))

    # Construct the problem.
    if not M_is_empty:
        X_left = cvx.Variable((num_nonzero_M, 1))
    Y_left = cvx.Variable((num_nonzero_rep, 1))
    Z_left = cvx.Variable((P.shape[0], tau_max*n))

    obj_left = - C_coeff*Y_left + cvx.sum(cvx.multiply(P, Z_left))
    if not M_is_empty:
        obj_left += W_coeff*X_left

    if M_is_empty:
        cons_left = [R_left == row_sum_matrix_Y*Y_left,
                      cvx.sum(Z_left, axis=0, keepdims=True).T == col_sum_matrix_Y * Y_left + Ru_left,
                      0 <= Y_left, 0 <= Z_left, Z_left <= PLen]
    else:
        cons_left = [R_left == row_sum_matrix_X*X_left + row_sum_matrix_Y*Y_left,
                      cvx.sum(Z_left, axis=0, keepdims=True).T == col_sum_matrix_X * X_left + col_sum_matrix_Y * Y_left + Ru_left,
                      0 <= X_left, X_left <= M_coeff, 0 <= Y_left, 0 <= Z_left, Z_left <= PLen]

    prob_left = cvx.Problem(cvx.Maximize(obj_left), cons_left)

    # Solve
    if solver == 'ECOS_BB':

t : pd.timestamp
            Timestamp for the optimization.
        """

        if t is None:
            t = pd.datetime.today()

        value = sum(portfolio)
        w = portfolio / value
        z = cvx.Variable(w.size)  # TODO pass index
        wplus = w.values + z

        if isinstance(self.return_forecast, BaseReturnsModel):
            alpha_term = self.return_forecast.weight_expr(t, wplus)
        else:
            alpha_term = cvx.sum(cvx.multiply(
                time_locator(self.return_forecast, t, as_numpy=True), wplus))

        assert(alpha_term.is_concave())

        costs, constraints = [], []

        for cost in self.costs:
            cost_expr, const_expr = cost.weight_expr(t, wplus, z, value)
            costs.append(cost_expr)
            constraints += const_expr

        constraints += [item for item in (con.weight_expr(t, wplus, z, value)
                                          for con in self.constraints)]

        for el in costs:
            assert (el.is_convex())

def _estimate(self, t, w_plus, z, value):
        """Estimate holding costs.

        Args:
          t: time of estimate
          wplus: holdings
          tau: time to estimate (default=t)
        """
        try:
            w_plus = w_plus[w_plus.index != self.cash_key]
            w_plus = w_plus.values
        except AttributeError:
            w_plus = w_plus[:-1]  # TODO fix when cvxpy pandas ready

        try:
            self.expression = cvx.multiply(
                time_locator(self.borrow_costs, t), cvx.neg(w_plus))
        except TypeError:
            self.expression = cvx.multiply(time_locator(
                self.borrow_costs, t).values, cvx.neg(w_plus))
        try:
            self.expression -= cvx.multiply(
                time_locator(self.dividends, t), w_plus)
        except TypeError:
            self.expression -= cvx.multiply(
                time_locator(self.dividends, t).values, w_plus)

        return cvx.sum(self.expression), []

no_trade = second_term.index[second_term.isnull()]
            second_term[no_trade] = 0
            constr += [z[second_term.index.get_loc(tick)] == 0
                       for tick in no_trade]

        try:
            self.expression = cvx.multiply(
                time_locator(self.half_spread, t), cvx.abs(z))
        except TypeError:
            self.expression = cvx.multiply(
                time_locator(self.half_spread, t).values, cvx.abs(z))
        try:
            self.expression += cvx.multiply(second_term,
                                            cvx.abs(z) ** self.power)
        except TypeError:
            self.expression += cvx.multiply(
                second_term.values, cvx.abs(z) ** self.power)

        return cvx.sum(self.expression), constr

def _estimate(self, t, wplus, z, value):
        F = locator(self.exposures, t)
        f = (wplus.T * F.T).T
        Sigma_F = locator(self.factor_Sigma, t)
        D = locator(self.idiosync, t)
        self.expression = cvx.sum_squares(
            cvx.multiply(np.sqrt(D), wplus)) + \
            cvx.quad_form(f, Sigma_F) + \
            self.epsilon * (cvx.abs(f).T * np.sqrt(np.diag(Sigma_F)))**2

        return self.expression

The objective is to maximize the total utility, sum(ui(Ri),i=1..n)
where ui: R → R is the utility function associated with the ith receiver.
  '''
  # Input parameters: alpha and beta are constants from R_i equation
  n=len(a_val)
  if n!=len(b_val):
    print('alpha and beta vectors must have same length!')
    return 'failed',np.nan,np.nan,np.nan
  P=cvx.Variable(n)
  W=cvx.Variable(n)
  alpha=cvx.Parameter(n)
  beta =cvx.Parameter(n)
  alpha.value=np.array(a_val)
  beta.value =np.array(b_val)
  # This function will be used as the objective so must be DCP; i.e. element-wise multiplication must occur inside kl_div, not outside otherwise the solver does not know if it is DCP...
  R=cvx.kl_div(cvx.multiply(alpha, W),
               cvx.multiply(alpha, W + cvx.multiply(beta, P))) - \
    cvx.multiply(alpha, cvx.multiply(beta, P))
  objective=cvx.Minimize(cvx.sum(R))
  constraints=[P>=0.0,
               W>=0.0,
               cvx.sum(P)-P_tot==0.0,
               cvx.sum(W)-W_tot==0.0]
  prob=cvx.Problem(objective, constraints)
  prob.solve()
  return prob.status,-prob.value,P.value,W.value