How can I simulate and display data that satisfies several constrains?

Hi there, I confess I am quite new to MatLab so the following might sound very ammateur... I need to simulate uniformly data (two values: x and y) that obey a constrain/condition. I was able to use the function rand but I cannot get an easy solution to my problem.... Here it goes...Suppose I have the following 28 constrains:

1) x+2y=100 2) x+3y=100 3) 2x+y=100 4) 2x+3y=100 5) 3x+y=100 6) x+y=100 7) 3x+2y=100 8) x+2y=40 9) x+3y=40 10) 2x+y=40 11) 2x+3y=40 12) 3x+y=40 13) x+y=40 14) 3x+2y=40 15) x+2y=80 16) x+3y=80 17) 2x+y=80 18) 2x+3y=80 19) 3x+y=80 20) x+y=80 21) 3x+2y=80 22) x+2y=60 23) x+3y=60 24) 2x+y=60 25) 2x+3y=60 26) 3x+y=60 27) x+y=60 28) 3x+2y=60

The above constrains can be generalized in the form: A.x+B.y=M . Now I need to generate 100 simulated values (i.e. 100 pairs of x and y) for EACH of these constrains. Using rand and previous solution from Matt Tearle I was able to do the following in MatLab:

function [y,x]=simula(B,A,M,n)
  y=M/B*rand(n,1);
  x=M/A -(B/A)*y;
  solution=[y,x];
  save solution solution
end

The above only gives me 100 pairs (n=100) of x and y for ONE constrain, so I need to do it for all of them and I need the solution in a user friendly layout. This means I need it in a matrix/table with the following 3 columns (and 2800 rows): i) column 1: the simulation number, which will be the same for 28 rows ( i.e. the 28 ordered constrains); ii) column 2: values for x (ordered according to the respective constrains); iii) column 3: values for y (ordered according to the respective constrains).

Many many thanks in advance. Kind regards, Paulo

[Walter Roberson]

Are the values generated required to solve all of the constraints simultaneously? Or you want 100 that satisfy #1, andother 100 that satisfy #2, and so on?

[Paulo]

Hi Walter, thanks for the question. The numbers need to solve the constrains NOT simultaneously. So, 100 numbers to solve #1, another 100 to solve #2, and so on. All the best

Answer 1

[Cedric]

If you are new to MATLAB and not looking for efficiency or elegance, you can use the following code that is easy to understand:

 a = [1 -3  7 -4] ;                    % Coefficients for x.
 b = [4  2 -1  3] ;                    % Coefficients for y.
 m = [0  2 -1  7] ;                    % Constant terms.

 nc = numel(a) ;                       % Number of "constraints".
 n  = 2 ;                              % Distribution size for x, y.

 results = zeros(nc*n, 3) ;            % Prealloc memory for results.

 for k = 1 : nc
    y = m(k)/b(k) .* rand(n, 1) ;      % Unif. distrib. for y.
    x = (m(k)-b(k).*y) ./ a(k) ;       % Distribution of x for given y.
    c = k .* ones(n, 1) ;              % Vector of current constr. ID.
    r = (k-1)*n + [1:n] ;              % Index range for current constr.
    results(r,:) = [c, x, y] ;         % Update block for current constr.
 end

In your case n=100 and a, b, m have to be adapted, but I kept them small here so you can check that it works. If the y distribution must be the same for all "constraints", you can just generate a vector of random numbers before the for loop statement and use it to define y, i.e.

 ...
 r = rand(n, 1) ;
 for k = 1 : nc
    y = m(k)/b(k) .* r ;
    ...

Here is a function with the same args order as in your question:

 function results = simula(b, a, m, n)
    % ...

    nc = numel(a) ;                    % Number of "constraints".
    results = zeros(nc*n, 3) ;         % Prealloc memory for results.

    for k = 1 : nc
       y = m(k)/b(k) .* rand(n, 1) ;   % Unif. distrib. for y.
       x = (m(k)-b(k).*y) ./ a(k) ;    % Distribution of x for given y.
       c = k .* ones(n, 1) ;           % Vector of current constr. ID.
       r = (k-1)*n + [1:n] ;           % Index range for current constr.
       results(r,:) = [c, x, y] ;      % Update block for current constr.
    end
 end

[EDITED @ 9:52EST] Finally, here is a slightly more efficient, yet not much more complicated version of it:

 C = ones(n, nc) ;
 X = zeros(n, nc) ;
 Y = rand(n, nc) ;
 for k = 1 : nc
    C(:,k) = k .* C(:,k) ;
    Y(:,k) = m(k)/b(k) .* Y(:,k) ;
    X(:,k) = (m(k)-b(k).*Y(:,k)) ./ a(k) ;
 end
 results = [C(:), X(:), Y(:)] ;

Note that C could be computed outside of the loop C=ones(n,nc)*diag(1:nc), but this is not as fast as computing it in the loop (in this context specifically).

Best regards,

Cedric

[Paulo]

Hi Cedric, MANY MANY THANKS for this! It's really nice and works really well. I just have a small request (sorry!): how can I have the solution grouped in sets of 28 constrains? This means for 100 simulations I would have 3 columns (and 2800 rows): i) column 1: sets of 28 rows each of which would have the 28 constrains); ii) column 2: the associated values for x (ordered according to the respective constraint); iii) column 3: the associated values for y (ordered according to the respective constraint). So if I run 2 simulations I would get 2 sets of 28 rows in 3 columns: column 1: contraint number 1 til 28 (twice); column 2: first 28 numbers are the associated values for x for simulation #1 and the following 28 numbers would be the values for x for simulation #2; and column 3: first 28 numbers are the associated values for y for simulation #1 and the following 28 numbers would be the values for y for simulation #2. Is this doable? Many thanks in advance. All the best Paulo

[Cedric]

Hi Paulo, you're most welcome. Just a question then, do you want the 28 constraints to be evaluated with the same random number or different ones? I.e. do you want this 100 times:

 r = rand(1) ;
 y(1) = m(1)/b(1) * r ;  x(1) = (m(1)-b(1)*y(1))/a(1) ;
 y(2) = m(2)/b(2) * r ;  x(2) = (m(2)-b(2)*y(2))/a(2) ;
 ...

Or this 100 times:

 r = rand(1) ;  y(1) = m(1)/b(1) * r ;  x(1) = (m(1)-b(1)*y(1))/a(1) ;
 r = rand(1) ;  y(2) = m(2)/b(2) * r ;  x(2) = (m(2)-b(2)*y(2))/a(2) ;
 ...

In the first case we pick a single random number and we evaluate the 28 constr. with it (same for all constr.), and we repeat that 100 times (so you have 100 experiments based on different random numbers, but within each experiment all constr. are evaluated based on the same random number). In the second case, each constraint is evaluated with a different random number.