converte the anonymous constraint function using fcn2optimexpr

N=20

N = 20

x=optimvar('x',N,'LowerBound',-1,'UpperBound',1);
y=optimvar('y',N,'LowerBound',-1,'UpperBound',1);
z=optimvar('z',N,'LowerBound',-2,'UpperBound',0);
elecprob=optimproblem;
cons1=x.^2+y.^2+(z+1).^2;
elecprob.Constraints.cons1=cons1<=1;

% elecprob.Constraints.plane1 = z <= -x-y;
% elecprob.Constraints.plane2 = z <= -x+y;
% elecprob.Constraints.plane3 = z <= x-y;
% elecprob.Constraints.plane4 = z <= x+y;

%%   if use anonymous constraint function ra[lace the code upside,
%%   the results are different,and there is something wrong with the 
%%   anonymous constraint function  and fcn2optimexpr
con2=@(x,y,z)(z+abs(x)+abs(y));
cons2=fcn2optimexpr(con2,x,y,z);
elecprob.Constraints.cons2=cons2<=0;
expr=optimexpr(1);
for ii=1:(N-1)
    for jj=(ii+1):N
       exp(ii)=(x(ii)-x(jj))^2+(y(ii)-y(jj))^2+(z(ii)-z(jj))^2;
       expr=expr+(exp(ii))^(-1/2);
    end
end

elecprob.Objective=expr;

rng default % For reproducibility
x0 = randn(N,3);
for ii=1:N
    x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;
    x0(ii,3) = x0(ii,3) - 1;
end
init.x = x0(:,1);
init.y = x0(:,2);
init.z = x0(:,3);
[sol,fval,flag]=solve(elecprob,init)

Solving problem using fmincon.

Feasible point with lower objective function value found.


Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3.000000e+03.

sol = struct with fields:
    x: [20×1 double]
    y: [20×1 double]
    z: [20×1 double]

fval = 166.2480

flag = 
    SolverLimitExceeded

If use anonymous constraint function replace the code in the % noted , the fval results 166.2480 are different from right result :163.0099, and there is something wrong with the anonymous constraint function and fcn2optimexpr.

[Torsten]

You shouldn't name a variable "exp" (see below) because exp is reserved for the exponential function in MATLAB. But the optimal value of the objective remains the same.

Are you sure about the optimal value being 163.0099 instead of 166.2480 ? Could you set the optimal x,y and z values as initial values and see what happens ?

N=20;
x=optimvar('x',N,'LowerBound',-1,'UpperBound',1);
y=optimvar('y',N,'LowerBound',-1,'UpperBound',1);
z=optimvar('z',N,'LowerBound',-2,'UpperBound',0);
elecprob=optimproblem;
cons1=x.^2+y.^2+(z+1).^2;
elecprob.Constraints.cons1=cons1<=1;

% elecprob.Constraints.plane1 = z <= -x-y;
% elecprob.Constraints.plane2 = z <= -x+y;
% elecprob.Constraints.plane3 = z <= x-y;
% elecprob.Constraints.plane4 = z <= x+y;

%%   if use anonymous constraint function ra[lace the code upside,
%%   the results are different,and there is something wrong with the 
%%   anonymous constraint function  and fcn2optimexpr
con2=@(x,y,z)(z+abs(x)+abs(y));
cons2=fcn2optimexpr(con2,x,y,z);
elecprob.Constraints.cons2=cons2<=0;
expr=optimexpr;
for ii=1:(N-1)
    for jj=(ii+1):N
       val=(x(ii)-x(jj))^2+(y(ii)-y(jj))^2+(z(ii)-z(jj))^2;
       expr=expr+val^(-1/2);
    end
end

elecprob.Objective=expr;

rng default % For reproducibility
x0 = randn(N,3);
for ii=1:N
    x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;
    x0(ii,3) = x0(ii,3) - 1;
end
init.x = x0(:,1);
init.y = x0(:,2);
init.z = x0(:,3);

opts = optimoptions(elecprob)
opts.MaxFunctionEvaluations = 10000
[sol,fval,flag]=solve(elecprob,init,Options=opts)
sol.x
sol.y
sol.z

Answer 1

[Matt J]

N=20;
x=optimvar('x',N,'LowerBound',-1,'UpperBound',1);
y=optimvar('y',N,'LowerBound',-1,'UpperBound',1);
z=optimvar('z',N,'LowerBound',-2,'UpperBound',0);
elecprob=optimproblem;
cons1=x.^2+y.^2+(z+1).^2;
elecprob.Constraints.cons1=cons1<=1;

%%   if use anonymous constraint function ra[lace the code upside,
%%   the results are different,and there is something wrong with the 
%%   anonymous constraint function  and fcn2optimexpr
con2=@(x,y,z)(z+abs(x)+abs(y));
cons2=fcn2optimexpr(con2,x,y,z);
elecprob.Constraints.cons2=cons2<=0;

[I,J]=find(triu(ones(N),1));

elecprob.Objective=sum(1./ sqrt( (x(I)-x(J)).^2 +   (y(I)-y(J)).^2   +  (z(I)-z(J)).^2)    );

rng default % For reproducibility
x0 = randn(N,3);
for ii=1:N
    x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;
    x0(ii,3) = x0(ii,3) - 1;
end
init.x = x0(:,1);
init.y = x0(:,2);
init.z = x0(:,3);

opts = optimoptions(elecprob);
opts.MaxIterations = 50000;
opts.OptimalityTolerance=1e-16;
opts.StepTolerance=1e-16;
opts.FunctionTolerance=1e-16;
opts.Algorithm='active-set';
[sol,fval,flag]=solve(elecprob,init,Options=opts);

Solving problem using fmincon.

Feasible point with lower objective function value found.


Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 6.000000e+03.

fval

fval = 162.9497

[Torsten]

@Matt J

You think the reason for the difficulty of the solver to converge is the non-differentiability of con2 ?

[jin yong]

how did you analysis this problem and dicided to set those options?

[Matt J]

@Torsten it probably does not help.

@jin yong Trial and error, mostly.

Answer 2

[jin yong]

the optimal value fval will be 163.0099.

replace the con2=@(x,y,z)(z+abs(x)+abs(y))

cons2=fcn2optimexpr(con2,x,y,z)

elecprob.Constraints.cons2=cons2<=0

with :elecprob.Constraints.plane1 = z <= -x-y;

elecprob.Constraints.plane2 = z <= -x+y;

elecprob.Constraints.plane3 = z <= x-y;

elecprob.Constraints.plane4 = z <= x+y;

the result will be: fval=163.0099

clear
N=20;
x=optimvar('x',N,'LowerBound',-1,'UpperBound',1);
y=optimvar('y',N,'LowerBound',-1,'UpperBound',1);
z=optimvar('z',N,'LowerBound',-2,'UpperBound',0);
elecprob=optimproblem;
cons1=x.^2+y.^2+(z+1).^2;
elecprob.Constraints.cons1=cons1<=1;
elecprob.Constraints.plane1 = z <= -x-y;
elecprob.Constraints.plane2 = z <= -x+y;
elecprob.Constraints.plane3 = z <= x-y;
elecprob.Constraints.plane4 = z <= x+y;
%%   if use anonymous constraint function ra[lace the code upside,
%%   the results are different,and there is something wrong with the 
%%   anonymous constraint function  and fcn2optimexpr
% con2=@(x,y,z)(z+abs(x)+abs(y));
% cons2=fcn2optimexpr(con2,x,y,z);
% elecprob.Constraints.cons2=cons2<=0;
expr=optimexpr;
for ii=1:(N-1)
    for jj=(ii+1):N
       val=(x(ii)-x(jj))^2+(y(ii)-y(jj))^2+(z(ii)-z(jj))^2;
       expr=expr+val^(-1/2);
    end
end
elecprob.Objective=expr;
rng default % For reproducibility
x0 = randn(N,3);
for ii=1:N
    x0(ii,:) = x0(ii,:)/norm(x0(ii,:))/2;
    x0(ii,3) = x0(ii,3) - 1;
end
init.x = x0(:,1);
init.y = x0(:,2);
init.z = x0(:,3);
opts = optimoptions(elecprob)

opts = 
  fmincon options:

   Options used by current Algorithm ('interior-point'):
   (Other available algorithms: 'active-set', 'sqp', 'sqp-legacy', 'trust-region-reflective')

   Set properties:
     No options set.

   Default properties:
                    Algorithm: 'interior-point'
           BarrierParamUpdate: 'monotone'
               CheckGradients: 0
          ConstraintTolerance: 1.0000e-06
                      Display: 'final'
        EnableFeasibilityMode: 0
     FiniteDifferenceStepSize: 'sqrt(eps)'
         FiniteDifferenceType: 'forward'
         HessianApproximation: 'bfgs'
                   HessianFcn: []
           HessianMultiplyFcn: []
                  HonorBounds: 1
       MaxFunctionEvaluations: 3000
                MaxIterations: 1000
               ObjectiveLimit: -1.0000e+20
          OptimalityTolerance: 1.0000e-06
                    OutputFcn: []
                      PlotFcn: []
                 ScaleProblem: 0
    SpecifyConstraintGradient: 0
     SpecifyObjectiveGradient: 0
                StepTolerance: 1.0000e-10
          SubproblemAlgorithm: 'factorization'
                     TypicalX: 'ones(numberOfVariables,1)'
                  UseParallel: 0

   Show options not used by current Algorithm ('interior-point')

opts.MaxFunctionEvaluations = 10000

opts = 
  fmincon options:

   Options used by current Algorithm ('interior-point'):
   (Other available algorithms: 'active-set', 'sqp', 'sqp-legacy', 'trust-region-reflective')

   Set properties:
       MaxFunctionEvaluations: 10000

   Default properties:
                    Algorithm: 'interior-point'
           BarrierParamUpdate: 'monotone'
               CheckGradients: 0
          ConstraintTolerance: 1.0000e-06
                      Display: 'final'
        EnableFeasibilityMode: 0
     FiniteDifferenceStepSize: 'sqrt(eps)'
         FiniteDifferenceType: 'forward'
         HessianApproximation: 'bfgs'
                   HessianFcn: []
           HessianMultiplyFcn: []
                  HonorBounds: 1
                MaxIterations: 1000
               ObjectiveLimit: -1.0000e+20
          OptimalityTolerance: 1.0000e-06
                    OutputFcn: []
                      PlotFcn: []
                 ScaleProblem: 0
    SpecifyConstraintGradient: 0
     SpecifyObjectiveGradient: 0
                StepTolerance: 1.0000e-10
          SubproblemAlgorithm: 'factorization'
                     TypicalX: 'ones(numberOfVariables,1)'
                  UseParallel: 0

   Show options not used by current Algorithm ('interior-point')

[sol,fval,flag]=solve(elecprob,init,Options=opts)

Solving problem using fmincon.

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.

sol = struct with fields:
    x: [20×1 double]
    y: [20×1 double]
    z: [20×1 double]

fval = 163.0099

flag = 
    OptimalSolution

sol.x

ans = 20×1

-0.5374 0.0000 -0.4736 -0.2989 0.5374 -1.0000 -0.7883 0.0000 0.7883 0.4736

sol.y

ans = 20×1

0.0305 -0.4877 0.0000 0.9271 0.0305 0.0000 0.5278 0.0000 0.5278 0.0000

sol.z

ans = 20×1

-1.8428 -1.8730 -0.4736 -1.2260 -1.8428 -1.0000 -1.3162 -0.0000 -1.3162 -0.4736