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
