Generate set of feasible points for a function?

Suppose you have a function (equality or inequality). How can you generate the set of feasible points from it?

For example: if you have 5x1 + 5x2 + 3x3 <= 8, the set of feasible points is {(0,0,0), (1,0,0), (0,1,0), etc.}. How can you use MATLAB to generate this set of points?

Answer 1

[Roger Stafford]

When asking for a set of all feasible points with respect to some condition to be generated by a computer you must ensure that there are only finitely many such points. In your example, there are infinitely many such points in the R^3 continuum. Even if x1, x2, and x3 are restricted to integer values, there are still infinitely many. If you add the condition that they each be an integer, say, in the range [0,10], then finally you have a finite number. A brute force method of generation in this last case would be:

 [X1,X2,X3] = ndgrid(0:10);
 X1 = X1(:); X2 = X2(:); X3 = X3(:);
 T = 5*X1+5*X2+3*X3<=8;
 A = [X1(T),X2(T),X3(T)];

It is impossible to make any general statement about an efficient method of generation. It would all depend on the nature of the particular condition imposed.

[Mike Vukovich]

Actually, that's indeed the case - I think that x1, x2 and x3 are restricted to positive integers. So in this case, would this be the way to generate the set of feasible points?

[Roger Stafford]

As I said, this would be a brute force method, but it should do the job. (Actually we could easily change the range to [0,8] since 9 and 10 are not possible values for any of the X's. It would still have to run through 9^3 = 729 possibilities.) It is possible that a for-loop method, properly done, could outperform the above, but the time required is rather small in any case. If you had a much larger value than 8 in your above condition, the number of trials by the brute force technique would increase greatly and a for-loop approach would surely be superior.

[Roger Stafford]

The for-loop method I mentioned might look like this:

 A = [];
 for x1 = 0:8/5
   for x2 = 0:(8-5*x1)/5
     t = (0:(8-5*x1-5*x2)/3)';
     A = [A;repmat([x1,x2],length(t),1),t];
   end
 end

As it stands this would take only three passes through the inner loop. With a larger number than 8 it would require substantially more but would likely be more efficient than the brute force method.

[Mike Vukovich]

Thanks! The brute-force method works well. However, the program I'm writing should be able to find the set of feasible points for any number of coordinates and constraints, not just 3; is there any way to do this?

[Roger Stafford]

It is easy to generalize the "brute force" technique to more than three variables and constraints, but not so easy to generalize methods that go faster. The latter would depend very much on the nature of those constraints. If there are n integer variables, x1, x2, x3, ...,xn and a set of m equalities (or in inequalities), eq1, eq2, ..., eqm, do this:

 [X1,X2,...,Xn] = ndgrid(a:b);
 X1 = X1(:); X2 = X2(:); ... ; Xn = Xn(:);
 T = eq1 & eq2 & ... & eqm; % <-- Each depends on X1, X2, ..., Xn
 A = [X1(T),X2(T), ...,Xn(T)];

1. The triple-dot ellipsis here refers to this text and is not the one used in matlab.

2. The quantities 'a' and 'b' are to be chosen as the least and greatest values that are possible for any X variable. If the variables have differing ranges, then 'ndgrid' should contain all these distinct ranges as arguments.

One note of warning. As the number of variables and their ranges increases, the total number of cases to be tested increases very rapidly and soon this brute force method becomes impractical. For example if you have 10 variables and each needs to be tested for 25 different possible integers, the total to be tested with this crude procedure would be 25^10 or nearly 10^14 = hundred million million. You would very likely have to devise some other more clever method.

[Mike Vukovich]

Do you know how I can generate the appropriate number of vectors depending on n? I tried this code, but regardless of what n is, only three vectors are generated.

Answer 2

[Azzi Abdelmalek]

a=arrangement([1 0],3);
co1=[5 5 3];
co2=8
b=sum(bsxfun(@times,a,co1),2)-co2;
out=a(b<=0,:)

% Before runing the above code, save the function arrangement.m:

function y=arrangement(v,n)
 m=length(v);
 y=zeros(m^n,n);
 for k = 1:n
    y(:,k) = repmat(reshape(repmat(v,m^(n-k),1),m*m^(n-k),1),m^(k-1),1);
 end

[Mike Vukovich]

Thanks, but how does it work? Also, not all the feasible points contain only 1 and 0, so would this function work for that?

[Azzi Abdelmalek]

Have you tried the codes