Minimizing 3x + x*e^(1/x) such that 1<= x <=10 using cvxpy

Question

I am trying to solve the following problem using CVXPY

import cvxpy as cp
x = cp.Variable(nonneg=True)
from cvxpy.atoms.elementwise.power import power
objective = cp.Minimize(3*x + x*cp.exp(cp.inv_pos(x)))
constraints = [x>=1,x<=10]
problem = cp.Problem(objective, constraints)
problem.solve()
print(x.value)
print(problem.value)

It gives me the following error:

DCPError: Problem does not follow DCP rules. Specifically:
The objective is not DCP. Its following subexpressions are not:
var91 @ exp(power(var91, -1.0))

I think this objective is convex, but there in an error

Answer 1

cvxpy

As @Michal Adamaszek mentioned, you can add another variable t saisfying the constraint cp.rel_entr(x, t) <= -1 i.e. t>=x*exp(1/x) (relaxation), then minimize 3*x+t

import cvxpy as cp

# Define the variables
x = cp.Variable(nonneg=True)
t = cp.Variable(nonneg=True)

# Define the objective function
objective = cp.Minimize(3*x + t)

# Define the constraints
constraints = [
    x >= 1,
    x <= 10,
    cp.rel_entr(x, t) <= -1  # Exponential cone constraint
]

# Define and solve the problem
problem = cp.Problem(objective, constraints)
problem.solve()

# Print the results
print("Optimal value of x:", x.value)
print("Optimal value of the objective function:", problem.value)

Optimal value of x: 0.9999999906618644
Optimal value of the objective function: 5.718281734146801

MATLAB, with CVX

clear all;close all;clc;

% Make sure to have CVX installed and set up

% Start a new CVX model
cvx_begin

    % Define the variables
    variable x nonnegative
    variable t nonnegative

    % Define the objective function
    minimize(3*x + t)

    % Define the constraints
    subject to
        x >= 1
        x <= 10
        rel_entr(x, t) <= -1  % Exponential cone constraint

% Solve the problem
cvx_end

% Print the results
disp(['Optimal value of x: ', num2str(x)]);
disp(['Optimal value of the objective function: ', num2str(cvx_optval)]);

Successive approximation method to be employed.
   For improved efficiency, SDPT3 is solving the dual problem.
   SDPT3 will be called several times to refine the solution.
   Original size: 8 variables, 3 equality constraints
   1 exponentials add 8 variables, 5 equality constraints
-----------------------------------------------------------------
 Cones  |             Errors              |
Mov/Act | Centering  Exp cone   Poly cone | Status
--------+---------------------------------+---------
  0/  0 | 0.000e+00  0.000e+00  0.000e+00 | Solved
-----------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +5.71828
 
Optimal value of x: 1
Optimal value of the objective function: 5.7183

MATLAB, with YALMIP

clear all;close all;clc;
% Define the decision variable
x = sdpvar(1);

% Define the objective function
objective = 3*x + x*exp(1/x);

% Define the constraints
constraints = [1 <= x <= 10];

% Setup the options for the solver, if necessary
% options = sdpsettings('solver', 'desired_solver_name');

% Solve the problem
sol = optimize(constraints, objective);

% Check if the solver found a solution and display it
if sol.problem == 0
    % Solution found
    disp('Optimal x:');
    disp(value(x));
else
    % Handle the case where no solution was found
    disp('No solution found.');
end

Your initial point x0 is not between bounds lb and ub; FMINCON
shifted x0 to strictly satisfy the bounds.

                                            First-order      Norm of
 Iter F-count            f(x)  Feasibility   optimality         step
    0       1    9.259209e+00    0.000e+00    1.907e+00
    1       2    5.733165e+00    0.000e+00    4.579e-01    9.850e-01
    2       3    5.834661e+00    0.000e+00    5.572e-02    3.322e-02
    3       4    5.721081e+00    0.000e+00    9.326e-04    3.723e-02
    4       5    5.718296e+00    0.000e+00    1.792e-04    9.279e-04
    5       6    5.718292e+00    0.000e+00    3.333e-06    1.330e-06
    6       7    5.718282e+00    0.000e+00    3.334e-08    3.300e-06

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

<stopping criteria details>
Optimal x:
    1.0000

gurobipy

from gurobipy import Model, GRB

# Create a new model
m = Model("optimizer")

# Create variable - lower bound is 1, upper bound is 10
x = m.addVar(lb=1, ub=10, name="x")

# Create a new auxiliary variable for the exponential term
exp_term = m.addVar(name="exp_term")

# Link the exponential variable with x using the exponential constraint
m.addGenConstrExp(x, exp_term, name="exp_constr")

# Set objective: minimize 3*x + x*exp(1/x)
# Since the exponential term is now represented by 'exp_term', we multiply it by x
m.setObjective(3*x + x*exp_term, GRB.MINIMIZE)

# Optimize the model
m.optimize()

# Print the optimal value of x
if m.status == GRB.Status.OPTIMAL:
    print('Optimal value of x: %.4g' % x.X)
else:
    print('Optimal solution was not found.')

Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (win64 - Windows 10.0 (19045.2))

CPU model: Intel(R) Core(TM) i7-10750H CPU @ 2.60GHz, instruction set [SSE2|AVX|AVX2]
Thread count: 6 physical cores, 12 logical processors, using up to 12 threads

Optimize a model with 0 rows, 2 columns and 0 nonzeros
Model fingerprint: 0xffe6a271
Model has 1 quadratic objective term
Model has 1 general constraint
Variable types: 2 continuous, 0 integer (0 binary)
Coefficient statistics:
  Matrix range     [0e+00, 0e+00]
  Objective range  [3e+00, 3e+00]
  QObjective range [2e+00, 2e+00]
  Bounds range     [1e+00, 1e+01]
  RHS range        [0e+00, 0e+00]
Presolve added 17 rows and 108 columns
Presolve time: 0.03s
Presolved: 19 rows, 111 columns, 1492 nonzeros
Presolved model has 1 bilinear constraint(s)

Solving non-convex MIQCP

Variable types: 104 continuous, 7 integer (1 binary)
Found heuristic solution: objective 5.7182818

Explored 0 nodes (0 simplex iterations) in 0.05 seconds (0.00 work units)
Thread count was 12 (of 12 available processors)

Solution count 1: 5.71828 

Optimal solution found (tolerance 1.00e-04)
Best objective 5.718281829554e+00, best bound 5.718281828456e+00, gap 0.0000%
Optimal value of x: 1

Mathematica

Plot[3 x + x*Exp[1/x], {x, 0.00001, 100}]
Minimize[{3 x + x*Exp[1/x], 1 <= x <= 10}, x, Reals]

{3 + E, {x -> 1}}