Particle swarm optimization solving an equation system

I want to use particle swarm optimization to solve 4 complex valued nonlinear equations just to compare the results I got from other optimizers. It is working but I'm getting a single value for each variable not in the form of "a+bi" or "a-bi" because of the norm I guess. When I remove the "norm" I am getting an error message. How can I fix the poroblem? Thanks for the help.

ObjFnc = @(x) [...]; % Equations
    
nVars = 4;
LB = [-1,-1,-1,-1];
UB = [1,1,1,1];

PS_Opts = optimoptions('particleswarm','HybridFcn',@fmincon); 
[x,fval,exitflag,output] = particleswarm(@(x)norm(ObjFnc(x)),nVars,LB,UB,PS_Opts)

Answer 1

[Dana]

As far as I know, you cannot directly optimize using complex variables as your optimizing variables, and of course your objective function must be a real scalar so you can't have a comlex variable there either.

All is not lost, however. Instead of having 4 complex variables x1, x2, x3, x4 as your optimizers, have 8 real variables a1, b1, a2, b2, a3, b3, a4, b4, and then in your objective function convert these to complex variables via x1=a1+b1*i, x2=a2+b2*i, etc.

[MByk]

Thanks for the answer Dana. I actually got the results as I wanted using "fsolve" and "lsqnonlin". I just want to try particle swarm optimization. I was thinking about using "real" and "imag" functions for each equation. If I got it right, you're suggesting the same solution. But this time I have to set the number of variable to 8, right?

[Dana]

If I understand correctly what it is you're doing, then I don't see the use for the real and imag functions.

As I understand it, you have a function ObjFnc that takes 4 complex variables x1,...,x4 as inputs, and outputs a complex vector with 4 elements in it. You seek the values of x1,...,x4 that makes the output of ObjFnc as close to the zero vector as possible, and a natural measure of this "closeness" is the norm of ObjFnc.

If that's correct, then you won't have any need for the real and imag functions. Those functions are useful for converting a single complex number into an equivalent pair of real numbers (i.e., the real and imaginary parts of the complex number). You're only going in the opposite direction here, though: you have real numbers (a1,b1,...,a4,b4) that you want to convert to complex numbers, which you can do simply via x1=a1+b1*i.

Actually, a more robust syntax is to use x1=a1+b1*1i, i.e., use 1i instead of i: while the default value of i is , i can in principle be overwritten with some other value at any time in a script or function (e.g., you can set i=3 if you want). On the other hand, 1i is always reserved to mean , so you don't have to worry about accidentally overwriting it with something.

In any case, yes, if you do as I suggested, you'd need to set the number of variables in the particle swarm to 8.

[MByk]

This is how our equation system looks like.

Sd1 = -0.063933409 - 0.18155469 * 1i;
...
Sd4 = -0.981188030 + 0.15023738 * 1i

ObjFnc = @(x) [
x(1)-Sd1-x(2).*((x(4).*(x(3)-x(2)+...;
-Sd2-x(4).*((x(2).*(x(3).*x(2)+x(3).*...;
x(1)-Sd3+x(2).*((x(2).*(x(4).^2.*x(2).*x(3).^2...;
-Sd4-x(4).*((x(4).*(x(2).*x(1)-1).^2)./(x(3).^2...];

[Dana]

Okay wait, so is x here complex or real? I had thought you were saying it was complex, and that's what my previous posts have been trying to show you how to deal with, but re-reading everything now I'm wondering if the x is actually real, and it's only the output of the function OjbFnc that's complex.

If that's the case, then I don't quite see what the problem you're having is. You should be able to do exactly what you wrote in your first post, i.e.,

[x,fval,exitflag,output] = particleswarm(@(x)norm(ObjFnc(x)),nVars,LB,UB,PS_Opts)

Yes, the output for the optimized function value fval would be a single real number, which is norm(ObjFnc(x)), but that makes sense: norm(ObjFnc(x)) is the function you're telling MATLAB to minimize. It also wouldn't make sense to just try to optimize the vector-valued ObjFnc without the norm: for two different values of x, MATLAB would have no idea how to compare ObjFnc(x) to see which is "smaller", since the concept of "smaller" isn't well-defined for vectors (and even if it were, you're not trying to minimize ObjFnc, you're trying to make it close to zero). This is why MATLAB threw an error when you removed the norm.

If all you want here is to know the actual values of the output of ObjFnc(x) at the optimum (rather than the norm of that value), all you need to do is call ObjFnc afterwards with the optimized value. So something like:

[xopt,fval,exitflag,output] = particleswarm(@(x)norm(ObjFnc(x)),nVars,LB,UB,PS_Opts);
ObjOpt = ObjFnc(xopt);

You can then inspect the elements of ObjOpt if you want. Is this what you're after?

[MByk]

Actually, we are trying to find complex valued xs.

ObjOpt = ObjFnc(xopt); works Dana. I am getting complex values now.

ObjOpt = ObjFnc(x);
ObjOpt =
   0.0036 + 0.1816i
   0.0027 + 0.3283i
   0.0004 + 0.1694i
   0.0014 + 0.3327i

Thank you so much for your assistance.

[Dana]

Just to make sure, if you're using complex-valued x's, then you need to use the original suggestion I gave you, where your optimization has 8 real variables as inputs, which you then use as the real and imaginary parts of your x's in your objective function. Otherwise particleswarm is only searching among real values of x (you can check this by noticing that your optimum values of x are real). Just because you're getting a complex output of ObjFnc, doesn't mean the inputs are complex: the constants Sd1,...,Sd4 that you're using are complex, which means the output of ObjFnc will be complex even if the x's are real.

[MByk]

Yes, you are right I missed that. I got what you mean but I don't know how to do it. Is there an example I can use?

[Dana]

First of all I think you'll find this easier to program if you make ObjFnc a proper function (i.e., in its own .m file), rather than an anonymous one as you have. Then that .m file could look something like:

function obj = ObjFnc(ab) 

% The input ab here would be a vector of length 8, with the first 4 elements 
% being the real parts of the x's, and the last 4 the complex parts.

Sd1 = -0.063933409 - 0.18155469 * 1i;
...
Sd4 = -0.981188030 + 0.15023738 * 1i

x = ab(1:4) + ab(5:8)*1i;

obj = x(1)-Sd1-x(2).*((x(4).*(x(3)-x(2)+...;
-Sd2-x(4).*((x(2).*(x(3).*x(2)+x(3).*...;
x(1)-Sd3+x(2).*((x(2).*(x(4).^2.*x(2).*x(3).^2...;
-Sd4-x(4).*((x(4).*(x(2).*x(1)-1).^2)./(x(3).^2...];

end

You can then optimize this with the particle swarm as

nVars = 8;

[abopt,fval,exitflag,output] = particleswarm(@(ab)norm(ObjFnc(ab)),nVars,LB,UB,PS_Opts);

The only remaining issue you'd need to address here is the constraints on the elements of ab. You originally had

LB = [-1,-1,-1,-1]; 
UB = [1,1,1,1];

but (a) that wouldn't make sense even if ObjFnc took complex inputs (is 0.9+0.2i < or > than 1?), and (b) you now need to put bounds on both the real and imaginary parts, so these vectors would need to each have 8 elements in them. If you don't actually have constraints, you can just omit these arguments altogether. If you do have constraints, you need to think carefully about how to deal with them, and depending on the form they take this may require a change of variables.

For example, if you want to constrain |xj| <= 1, then you're going to have to take a different approach, since this is a joint restriction on aj and bj, and you can't do that with particleswarm. In that case, rather than using the real and imaginary parts of the x's as the function inputs, you might be better off working in polar coordinates:

function obj = ObjFnc(thR) 

% The input thR here would be a vector of length 8, with th = thR(1:4), R = thR(5:8),
% and x(j) = R*exp(th*i) = R*(cos(th) + i*sin(th)).

Sd1 = -0.063933409 - 0.18155469 * 1i;
...
Sd4 = -0.981188030 + 0.15023738 * 1i

th = thR(1:4);
R = thR(5:8);
x = R*exp(th*1i);

obj = x(1)-Sd1-x(2).*((x(4).*(x(3)-x(2)+...;
-Sd2-x(4).*((x(2).*(x(3).*x(2)+x(3).*...;
x(1)-Sd3+x(2).*((x(2).*(x(4).^2.*x(2).*x(3).^2...;
-Sd4-x(4).*((x(4).*(x(2).*x(1)-1).^2)./(x(3).^2...];

end

In this case, you can bound the elements of R between -1 and 1, which will ensure that |x|<=1. Technically you could also bound the angle th between, say, and π, but this may create numerical issues near the boundaries. It might be better to just leave the angle unbounded (i.e., use bounds -Inf, Inf).

[MByk]

Wow. :) It's more complicated than I expected. Actually, I tried to modify anonymous function and failed. Transforming it to a proper function is a great idea. I am going to think about how to deal with the bounds. I have to ask my advisor. I greatly appreciate your help. I learned much from your detailed answers.

Answer 2

[Ceyhan TURKMEN]

Are the first 4 values real and the last 4 values imaginary part of the results?