Updating multiple variables randomly with a constraint in a while cycle

I have a pretty complex function that calculates a 1 by 6 vector (phi) from another 1 by 6 vector (x) like (phi)=function(x) I want now to fix phi and look for all the possible x that give me that phi with the costraint that sum(x)==1 and that for every element of x(i), 0<= x <= 1

I'm trying with a WHILE loop but I don't know how to randomly update the x variable preventing it to try the same x multiple times; my desired solution have 3 digits of accuracy like x=[0.0112 0.0997 0.157 0.227 ecc.] I know there could be a maximum of 2 solutions and I need to find them both eventually.

I don't necessarily need to use a WHILE loop If there is any MATLAB functionality that can solve this it would be great.

[gabriele provenza]

please ask me for more details if needed

[darova]

Please attach the code

[gabriele provenza]

Tc=[304.2,735.8,757.9,782.4,817.1,837.3,];                  %Kelvin

Pc=[7.376,1.454,1.315,1.236,1.269,1.046];      %MPa
w=[0.225,0.855,0.928,0.994,1.046,1.108];   %CO2 EEC14	EEC16	EEC18	EPA	 DHA
MW=[44,249,276,303,330,328.5];
P=11;
T=313;


alfaij=[0       0.07	-0.155	-0.208	-0.348	-0.036
        0.07    0       0       0       0       0
        -0.155  0       0       0       0       0
        -0.208  0       0       0       0       0
        -0.348  0       0       0       0       0
        -0.036  0       0       0       0       0];

betaij=[0       0       0.0007	0.0009	0.0013	0.0004
        0       0       0       0       0       0
        0.0007	0       0       0       0       0
        0.0009	0       0       0       0       0
        0.0013	0       0       0       0       0
        0.0004	0       0       0       0       0];

nij=[0          0       0.030	0.030	0.020	0.012
     0          0       0       0       0       0
     0.030      0       0       0       0       0
     0.030      0       0       0       0       0
     0.020      0       0       0       0       0
     0.012      0       0       0       0       0];
 
x=[0.001,0.2,0.2,0.2,0.2,0.19];
phix=PENG_FUGACITY_function(P,T,1,Tc,Pc,w,x,alfaij,betaij,nij);

y=[0.99,0.005,0.005,10^-15,10^-15,10^-15];
    
while abs(sum(abs(phiy-phix))>0.001


phiy=PENG_FUGACITY_function(P,T,1,Tc,Pc,w,y,alfaij,betaij,nij);

%here i should update y with a certain logic i'm looking for

[darova]

• I want now to fix phi and look for all the possible x that give me that phi with the costraint that sum(x)==1;

Maybe you need fsolve or fmincon. Looks like optimization problem

[gabriele provenza]

i'm trying to grasp how to use fmincon, hope it will work; unfortunately in this case there should be 2 solutions and I'd need a specific one of them (or both) I'll see what i can do.

Thanks

[John D'Errico]

Incredibly confusing question. However you cannot use fmincon to choose anything randomly. You cannot use fmincon to solve a discrete problem, which this somehow seems like it might be. And finally, you cannot get fmincon to return two solutions. And of course we don't see your code, but only how you call it.

You can use fmincon with a constraint that the sum of the unknowns is a contstant, thus 1 here. But the rest of what you are asking for is really confused.

It is not clear how it is that you are sure there are exactly two solutions either.

Are you trying somehow to perform some sort of ad hoc optimization without the use of fmincon? Is that what the idea of choosing values randomly is intended to do?

[gabriele provenza]

I know there are only 2 solutions because this modellization is about a two phases equilibrium, already studied by other persons, bibliography suggests there should be only 2 solutions.

Given that the accuracy needed is 3 digits and given the presence of the constraints it should be possible to compute every combination of 6 numbers such that the sum is one and (I updated the question since i forgot to mention it) every number can vary between 0 and 1. As i was saying the objetive was to try every combination; how could i do so? I tought there could be some random generator that doesn't produce duplicates. If there is any way to let a function calculate an optimum solution it would be good either.

I was trying to discretize in order to be able to perform the calculation for every case and check only the cases that give me a very low error; I never used such functions and tought I should have somehow to work with a while cycle, so my problem was about logic; with one variable varying between 0 and 1 I could go for something like

while y-5>0.00001

x=0;

y=f(x)

x=x+0.00001

if x>1 break

end

end

how can I update the x in such a way that it covers the whole domain of my constraints if I have a 1 by 6 array as x instead of a variable?

[Sindar]

While I don't know a simple way to iterate through the possibilities, counting them is straightforward. This is a classic "stars and bars" problem: you have 1000 units of 0.001 to split between the 6 elements. The number of possibilities is then (1000-1) choose (6-1), or about 8e12. This is very probably more than you can reasonably calculate.

Note: it's not clear what exactly you mean by "3 digits of accuracy," since your examples x elements have 4 digits after the decimal (3 significant figures?) and your last comment shows x with 5 digits. Still, it only gets (much much) worse if you have finer-than-0.001 precision

[gabriele provenza]

yes, i mean 3 significant figures, some number could reasonably be in the order of 10^-6

[Sindar]

In that case, the number of possible solutions is dramatically more, and significantly harder to either count or iterate over

[Sindar]

Therefore, my suggestion is to drop the discretization, and treat it as a continuous optimization problem. Once you have an answer, truncate it to your desired accuracy and check how poor the agreement is.

Answer 1

[Jeff Miller]

You might be able to do what you want with nested while loops, something like this. Pretty tedious though.

x = zeros(6,1);
step = 0.001;  % minimum step between adjacent x values
x1max = 1 - 5*step;  % x1 cannot be bigger than this because at least step is needed for all others
x(1) = 0;
nFound = 0;  % number of solutions found so far
while x(1)<=x1max && nFound<2
    x(1) = x(1) + step;
    x2max = 1 - x(1) - 4*step;
    x(2) = 0;
    while x(2)<=x2max && nFound<2
       x(2) = x(2) + step;
       x3max = 1 - x(1) - x(2) - 3*step;
       x(3) = 0;
       % add further nested while loops like the above for x(3), x(4) and x(5)
       % finally, deep in the midst of all these while loops, do this
                  x(6) = 1 - sum(x(1:5));
                  phi = func(x);
                  % check phi here.
                  % If phi is a solution you want, save it and increment nFound
    end
end

[Sindar]

While this technically would work, there are several issues for this particular problem:

• the number of iterations is ridiculous (I estimated 8e12 above)
• Gabriele eventually clarified that they mean 3 significant figures, not 3 places after the decimal
• while possible, I think it unlikely that there is a way to check phi such that only the two desired solutions pass. Instead, I expect that this is a normal convergence problem where the closer x gets to x_solution, the closer phi gets to phi_solution. So, depending on what tolerance you set, you could either miss x_solution ("nah, this phi is 1e-16 off from the real one") or find solutions that are not distinct (one step apart)
• since you are iterating through in a set order, it doesn't make much sense to use "while" over "for". It makes it harder to track progress, restart from where you left off, parallelize, etc..

Still, this is close to the optimal solution to the particular request given. But, it won't work. Conclusion: use a different algorithm