linear equation with restrictions
Leon galushko
i would carry out the following equation
(figures are as examples):
A1 * c1 + A2 * c2 = b
A = [2.30 4.56]
b = 57.4
the core problem is I'am looking for c1 and c2 which
can/should assume only integer values and are declared in an
array c = int32[(d e)]
where d = int32(1:1000) AND e = int32(1:1000)
Now the term to get computing started:
c = A\b
to solve the equation where c is restricted to an array of
given integers (1 to 10000)does not involve the restrictions
I had given to refer to by making computations and gives me
solving like this:
c =
0
12.5877
I'would be glad to have tips to solve this!
Bruno Luong
message <g5sgla$b7r$1...@fred.mathworks.com>...
> Hi,
> i would carry out the following equation
> (figures are as examples):
>
> A1 * c1 + A2 * c2 = b
> A = [2.30 4.56]
> b = 57.4
OK, in order to avoid work with float, I'll transform your
fisrt equation by multiplying it by 100
Ax100_1 * c1 + Ax100_2 * c2 = bx100 (eqt 1)
Ax100 = [230 456]
bx100 = 5740
Now search for greater common divisor of Ax1000 by MATLAB
function GCD:
[g D1 D2]=gcd(Ax100(1),Ax100(2))
% it is gcd(230,456) % g ==2, D1=-113, D2=57
Which mean
All the numbers generated by the lhs of (eqt 1) is
+ multiple of 2 (=g), and furthermore
+ g = 2 = -113*230 + 57*456 =
-113*Ax100(1) + 47*Ax100(2) (eqt 2).
Now all you need is to see whereas bx100 is divisible by g.
if mod(bx100 ,g)==0 % yes
C1 = D1 * bx100/g % -324310
C2 = D2 * bx100/g % +163590
else
error('No integer solution exists')
end
Verification:
C1*2.30+C2*4.56 % 57.4000
Bruno
Leon galushko
thanks indeed for your proposal to solve this!
I'm beginner of MatLab und therofore very thankfull for
initial support.
Now i have a question to this: g so far as i grasped is the
number of factors which are in a equation (C1, C2). Why is
'all the numbers generated by the ihs of (eqt1) must be
multiple of 2 (=g)'? And what is ihs?
Question2: c1 and c2 must be onle positiv, c1, c2 > 0.
How can this be secured?
Thanks!
Leon
"Bruno Luong" <b.l...@fogale.fr> wrote in message
<g5sjfa$fbg$1...@fred.mathworks.com>...
Bruno Luong
> Why is
> 'all the numbers generated by the ihs of (eqt1) must be
> multiple of 2 (=g)'? And what is ihs?
This is a well known theorem in algebra, it saids that all
the ring defined by {c1*a1 + c2*a2; c1, c2 in N} is an
"principal ideal" and the generator is the gcd of a1, a2.
> Question2: c1 and c2 must be onle positiv, c1, c2 > 0.
> How can this be secured?
In your case it cannot, or in other word there is no
solution that satisfies this constraint.
Here is why: If {d1,d2} is one couple integer solution (as
provided by gcd function) of (eqt2):
d1*a1 + d2*a2 = g, (eqt2)
then all solutions of (eqt 2) are (e1,e2) with
e1 = d1 + k*a2/g
e2 = d2 - k*a1/g
for k=..., -1, 0, 1, 2, ....
In your example, none of the solution has both (e1,e2)
positive. So (c1,c2) must have alternated sign.
Other example this might be possible, but not in your example.
Bruno
Leon galushko
Well, i might be bold and turn my qeustion more precisely,
indeed sorry for having failed this previosly!
The example with given figures means only to have made
problem clearly.
The assignment for MatLab would be to get only solvings
outputed which are positiv. If integer and positiv solutions
cant be found, then should be equally given a message like
you did in the code before: else
> error('No (positive)integer solution exists')
Probably there is a function for that or it can be done
through conditional looping, anyway you might have a right idea.
Leon
Bruno Luong
> then should be equally given a message like
> you did in the code before: else
> > error('No (positive)integer solution exists')
>
It's not that difficult, all we have to do is find a bracket
of an interval for 'k' so that both e1 and e2 are positives.
If the bracket contains at least one integer k, then there
is a solution (pick the k).
In the above example the bracket interval is
[0.495614035087719 0.495652173913044], and of course no
integer k belongs to it.
I let you figure out how to program this part in MATLAB.
Bruno
Leon galushko
Well done!
my really last question to:
in case there are more variables than 2 in the equation
(a1 * b2 + a2 * b2 + a3 * b3 = c) the solving with refering
to GCD - Function of MATLAB falls down becaus gcd operates
only with 2 elements?
Thanks!
Leon
Bruno Luong
No big deal, first find the gcd (g12) of the first two
coefs, then find the gcd (g123) of the g12 and a3, etc...
Bruno
Leon galushko
if i find gcd for a1 and a2 by this (refering to your code):
[g D1 D2]=gcd(Ax100(1),Ax100(2))
> % it is gcd(230,456) % g ==2, D1=-113, D2=57
[g ? D3]=gcd(Ax100(?),Ax100(3))
Bruno Luong
[g12 F1 F2]=gcd(Ax100(1),Ax100(2))
[g F12 F3]=gcd(g12,Ax100(3))
I little calculation is required to compute D1, D2, D3 from
F.
Bruno
Roger Stafford
> Well done!
> my really last question to:
> in case there are more variables than 2 in the equation
> (a1 * b2 + a2 * b2 + a3 * b3 = c) the solving with refering
> to GCD - Function of MATLAB falls down becaus gcd operates
> only with 2 elements?
> Thanks!
> Leon
Leon, I understand that you wish to solve the equation
a1*b1 + a2*b2 + a3*b3 = c
where a1, a2, a3, and c are known and without loss of generality can be
assumed to be positive integers, and where b1, b2, and b3 are unknown but
are restricted to also being positive integers.
If
[g12,-,-] = gcd(a1,a2)
is done, you already know how to find b1 and b2 which can make
a1*b1 + a2*b2
equal to any desired multiple of g12. If
[g123,-,-] = gcd(g12,a3)
is done, you can then find from this a set of b1, b2, and b3 values which can
make
a1*b1 + a2*b2 + a3*b3
equal to any desired multiple of g123. If c is not a multiple of g123, then the
original equation above has no solution. Otherwise b1, b2, and b3 solutions
can always be found, but one or two of them will typically be negative.
It is not hard to also show from the above reasoning that b3 can be adjusted
up or down by any desired multiple of g12/g123 if corresponding changes
are made to b1 and b2 that cause the left side of the equation to remain
unchanged. That is, if the multiple of g12 that a1*b1+a2*b2 is equal to is
changed in the opposite direction by the same multiple of a3/g123.
This should then be done in such a way that b3 assumes its minimum
possible positive value, and that will give c-a3*b3 its maximum possible
value. From what Bruno has told you, you already know how to then look for
possible solutions to the two variable problem
a1*b1 + a2*b2 = c3 = c-a3*b3
where b3 has been so determined. This will tell you if any valid positive
solutions exist and if so, will give you one of them.
If b3 were to be increased from this minimum value, then c-a3*b3 is
decreased and the number of possible positive b1 and b2 solutions can only
decrease. If you also have an upper bound for your b's, as in your first
article, the latter process of increasing b3 from its minimum might
conceivably be necessary to keep all b's within upper bounds simultaneously.
Roger Stafford
Leon galushko
wrote in message <g5ujb2$7c1$1...@fred.mathworks.com>...
thanks to your answer. I'm afraid i have but to clarify some
detail which Bruno almost tryed to satisfy.
(It's going about the same equation
(a1 * b2 + a2 * b2 + a3 * b3 = c) you understand rightly)
Say, if g1 and g2 through gcd(a1, a2) have been got
(output would be D1 and D2), then
to get g3 it's to state: gcd(g12,a3), where i understand
g12 to be Greatest Common Deviser for A1 and A1 (g12).
Where is have to be deduced g12 from?
Is it correct:
g12 = D1*A1 + D2*A2 ?
Thanks
Leon
Bruno Luong
Bruno
function [g varargout] = gcdn(varargin)
% function g = gcdn(a1, a2, ..., an); OR
%
% return g, Greatest common divisor of a1, ... an
%
% [g c1 c2 ... cn]=gcdn(a1, a2, ..., an)
% return also c1, ..., cn
% So that a1*c1 + ... an*cn = g
%
% Compact calling
% [g c]=gcdn(a1, a2, ..., an) or [g c]=gcdn(a)
% assumes a and c are array
%
if nargin<2
a = varargin{1};
if length(a)<2
error('gcdn requires two or more inputs');
end
a = reshape(a,1,[]);
else
% Put all numbers in array
a=cell2mat(varargin);
end
g=a(1);
c=zeros(size(a));
c(1)=1;
for k=2:length(a)
[g cg c(k)]= gcd(g, a(k));
c(1:k-1) = c(1:k-1)*cg;
end
if nargout>=2
switch (nargout-1)
case 1,
varargout{1}=c;
case length(c)
varargout=num2cell(c);
otherwise
error('number of output is incompatible with
input');
end
end
Bruno Luong
a = [75 180 35];
[g c]=gcdn(a)
sum(c.*a) % is g
Bruno
Roger Stafford
> Roger,
> thanks to your answer. I'm afraid i have but to clarify some
> detail which Bruno almost tryed to satisfy.
> (It's going about the same equation
> (a1 * b2 + a2 * b2 + a3 * b3 = c) you understand rightly)
> Say, if g1 and g2 through gcd(a1, a2) have been got
> (output would be D1 and D2), then
> to get g3 it's to state: gcd(g12,a3), where i understand
> g12 to be Greatest Common Deviser for A1 and A1 (g12).
>
> Where is have to be deduced g12 from?
> Is it correct:
>
> g12 = D1*A1 + D2*A2 ?
> Thanks
> Leon
I'll give a specific example, Leon, of what I was describing in the earlier
article to clarify things.
Suppose we have a1 = 770, a2 = 105, a3 = 60, and c = 2215, so we are
trying to solve
770*b1 + 105*b2 + 60*b3 = 2215.
First we do
[g12,d1,d2] = gcd(a1,a2)
which yields g12 = 35, d1 = 1, and d2 = -7. Next we do
[g123,e12,e3] = gcd(g12,a3)
which produces g123 = 5, e12 = -5, and e3 = 3. We notice that c is indeed a
multiple of g123: c = 2215 = 443*5, so we proceed. From above we know
that
e12*g12+e3*a3 = g123 [-5*35+3*60=5]
Therefore
(443*e12)*g12+(443*e3)*a3 = 443*g123, or
-2215*g12 + 1329*a3 = 2215
We can subtract any multiple of g12/g123=7 from 1329 if we add the same
multiple of a3/g123=12 to -2215, so we use the multiple 189 which gives
53*g12 + 6*a3 = 2215
This is the least positive value b3 could possibly have.
Hence we now wish to solve
a1*b1 + a2*b2 = c-6*a3 = 1855 = 53*g12
From above we know that
d1*a1 + d2*a2 = g12 [1*770-7*105=35]
so
(53*d1)*a1 + (53*d2)*a2 = 53*g12, or
53*a1 - 371*a2 = 1855
We can add any multiple of a1/g12=22 to -371 provided we subtract the
same multiple of a2/g12=3 from 53, so we choose 17 which gives
2*a1+3*a2 = 1855
Therefore we have a solution with all positive b's:
b1 = 2, b2 = 3, and b3 = 6.
Now you need to translate all the above reasoning into a general form that
can handle any situation.
I hope this example has been of assistance to you.
Roger Stafford
Bruno Luong
natural-integer linear equations with three or more
variables is a fun problem. Personally I would use a
recursive algorithm for it.
I'm not sure if there is a convenient way to parametrize the
set of infinite solutions. It would be nice if someone who
works with in number group/ring theory or
integer-programming can shed a light.
Bruno
Bruno Luong
of the integer equations:
given a: (1 x n), b : (1 x 1)
solve x integer, (1 x n) such that
a*x' = b
x >= 0
Example:
a = [770 105 60] % Like Roger's example
c = 7755
x =
0 3 124
0 7 117
0 11 110
0 15 103
0 19 96
0 23 89
0 27 82
0 31 75
0 35 68
0 39 61
0 43 54
0 47 47
0 51 40
0 55 33
0 59 26
0 63 19
0 67 12
0 71 5
3 1 89
3 5 82
3 9 75
3 13 68
3 17 61
3 21 54
3 25 47
3 29 40
3 33 33
3 37 26
3 41 19
3 45 12
3 49 5
6 3 47
6 7 40
6 11 33
6 15 26
6 19 19
6 23 12
6 27 5
9 1 12
9 5 5
Not 100% sure whereas my program is BUG free and would not
miss any solution.
Bruno
Bruno Luong
> I think finding all solutions (if they are finite) of this
> natural-integer linear equations with three or more
> variables is a fun problem.
No one else wants to tackle the integer problem "Help me to
spend all my money")? Roger?
I would like to see good ideas from you in solving this
problem (mine is rather brute-force).
I propose this a test example for us to compare notes:
a = [770 105 60 85 50]
b = 9435
Find *all* x (5x1), integer arrays such that
a*x= b and x>=0.
Bruno
Leon galushko
Bruno,
you said you made a program for outputting a set of all
possible combinations of solvings.
could you show your sourcecode of that program?
Leon
Bruno Luong
>
> could you show your sourcecode of that program?
>
Here we go, here is the main function, to be saved in one
m-file. The test program will follow.
function x = intlin(a, b)
% function x = intlin(a, b)
%
% PURPOSE: solving integer equation:
%
% x(1)*a(1) + x(2)*a(2) + ... + x(n)*a(n) = b
% x >=0
%
% Input: a (1 x n) integer (could be negative)
% b (1 x 1) integer (could be negative)
% Output: x (n x 1) integer such that
% a*x = b, x>=0
%
% Author: Bruno Luong
% Last update: 24/July/2008
% Default value, empty solution
x = zeros(0,length(a));
if ~isscalar(b)
fprintf('b must be scalar\n');
return
end
% cast to double, and reshape as long thin array
a = double(a(:));
b = double(b);
if ~all(mod(a,1)==0) || mod(b,1)~=0
fprintf('a and b must be integers\n');
return
end
% d*a' = g
[g d] = gcdn(a);
if mod(b,g)==0
an = a/g;
bn = b/g;
% (d*a' = b) OR (d*an' = bn)
d = bn*d;
%
% General gcd final solution would be
% x = d + k;
% with k such that: k*an' = 0, and
% k>=-d (<=> x>=0)
kmin = ceil(-d);
kmax = +inf(size(a));
% Find all k such that k.an = 0
k = allintlin0(an, kmin, kmax);
x = bsxfun(@plus, d, k);
else % mode(b,g)~=0
fprintf('WARNING: there is no solution\n');
end
end
function x = allintlin0(a, lower, upper)
%
% x, a, lower, upper are n-dimensional vector
% List all interger x such that
% x.a = 0
% lower <= x <= upper
%
a = a(:);
n = length(a);
lower = lower(:);
upper = upper(:);
% Adjust upper and lower bounds by LP
L = lower;
U = upper;
b = 0;
epsilon = 1e-6;
for k=1:n
cost = basis(k,n);
% Beware, his is Bruno's linprog, not MATLAB one in
optimization
% tool box, the result may be different.
sol = linprog([], [], a', b, cost, L, U);
if ~all(isinf(sol))
L(k) = max(L(k),sol(k)-epsilon);
end
cost = -basis(k,n);
sol = linprog([], [], a', b, cost, L, U);
if ~all(isinf(sol))
U(k) = min(U(k),sol(k)+epsilon);
end
end
L = ceil(L);
U = floor(U);
if all(~isinf(L)) && all(~isinf(U))
maxcount = Inf;
else % Limit the number of solutions that will be listed
% NOTE: set maxcount to finite value doesn't work as
efficienly,
% as expected because the recursive engine might spend
much of
% CPU time to look for unvalid solutions
% maxcount = 100;
fprintf('There is infinity of solutions\n');
x = NaN;
return
end
% Call the engine
count = 0;
x = ilinengine(a, L, U, b, count, maxcount);
if maxcount<inf
fprintf('WARNING: only %d solutions will be provided\n',
maxcount);
x = x(1:min(maxcount,end),:); % clipping
end
end
function v = basis(k,n) % generate a k-th basis vector of
dimension n
v = zeros(1,n);
v(k) = 1;
end
% Solver engine for integer x
% a*x = b
% lower <= x <= upper
% RESTRICTION:
% a must be primary array, i.e., they greatest common
divisor is one
function x = ilinengine(a, lower, upper, b, count, maxcount)
% default value, empty result
x = zeros(0,length(a));
if count>maxcount
return
end
% Trivial solution with one variable
if length(a)==1
if mod(b,a(1))==0
xtmp = b / a(1);
if (xtmp >= lower) && (xtmp <=upper)
x = xtmp;
end
end
return
end
% preliminary check where as the sum b is possible
% This check is to speed up (?), and deos not affect the result
as = bsxfun(@times,a,[lower upper]);
as = sum(sort(as,2),1);
if b<as(1) || b>as(2)
return
end
% Greatest common divisor for the tail
g = gcdn(a(2:end));
% find r1, the inverse of a(1) in g-modulo group
[uno r1] = gcd(a(1),g);
clear uno; % should be 1
r=r1*mod(b,g);
s = (b - a(1)*r)/g;
a(2:end) = a(2:end)/g; % the tail is primary among them
kmin = ceil((lower(1)-r)/g);
kmax = floor((upper(1)-r)/g);
% perform a basic step
function newcount = kstep(k)
% Recursive call
xk = ilinengine(a(2:end), lower(2:end),
upper(2:end), ...
s-k*a(1), count, maxcount);
nxk = size(xk,1);
x = [x; ...
(r+g*k)+zeros(nxk,1) xk]; % append the new solutions
newcount = count + nxk; % adjust the counter
end
if isinf(kmin) && isinf(kmax) % k is unbounded
for absk=0:inf
for k=unique([-absk absk])
if kstep(k)>maxcount
break
end
end
if count>maxcount
break
end
end
elseif isinf(kmin) % k has upper bound, but no lower bound
for k=kmax:-1:-inf
if kstep(k)>maxcount
break
end
end
else % k has lower bound
for k=kmin:kmax
if kstep(k)>maxcount
break
end
end
end
end
function [g varargout] = gcdn(varargin)
% function g = gcdn(a1, a2, ..., an);
%
% Return g, Greatest common divisor of a1, ... an
%
% [g c1 c2 ... cn]=gcdn(a1, a2, ..., an)
% return also c1, ..., cn
% So that a1*c1 + ... an*cn = g
%
% Compact calling form:
% [g c]=gcdn(a1, a2, ..., an) or [g c]=gcdn(a)
% assumes a and c are array
%
if nargin<2
a = varargin{1};
if length(a)<2
if a
g = abs(a);
if nargout>=2
varargout{1}=sign(a);
end
return
else
error('gcdn cannot compute for a = 0');
end
end
a = reshape(a,1,[]);
else
% Put all numbers in array
a=cell2mat(varargin);
end
g=a(1);
c=zeros(size(a));
c(1)=1;
for k=2:length(a)
[g cg c(k)]= gcd(g, a(k));
c(1:k-1) = c(1:k-1)*cg;
end
if nargout>=2
switch (nargout-1)
case 1,
varargout{1}=c;
case length(c)
varargout=num2cell(c);
otherwise
error('number of output is incompatible with
input');
end
end
end
Bruno Luong
a=[770 105 60 85 50];
b = 9435;
%
% solving
% c1*a1 + ... cn*an = b
% c >=0
c = intlin(a, b);
if isnan(c) % unbounded, infinity solutions
elseif ~isempty(c)
fprintf('\nSolutions c:\n\n');
disp(c);
fprintf('\nNumber of solutions = %d\n\n', size(c,1));
bverif = c*a';
bverifu = unique(bverif);
if length(bverifu)~=1 || bverifu~=b
fprintf('Wrong solution(s)\n');
disp(c(bverif~=b,:));
save test_inlin_debug.mat a cc b
end
else % mode(b,g)==0
fprintf('There is no solution\n');
end
Leon galushko
message <g6ic71$jqt$1...@fred.mathworks.com>...
Thank you for that!
I'am just trying to test it, now it comes to the following
error message:
??? Error using ==> linprog at 175
The number of rows in A must be the same as the length of b.
In your sourcecode I found a comment:
'% Beware, his is Bruno's linprog, not MATLAB one in
optimization
% tool box, the result may be different.
sol = linprog([], [], a', b, cost, L, U);'
Wouldn't run it with some Stadard- MatLab of mine?
Leon
Bruno Luong
>
> '% Beware, his is Bruno's linprog, not MATLAB one in
> optimization
> % tool box, the result may be different.
> sol = linprog([], [], a', b, cost, L, U);'
>
> Wouldn't run it with some Stadard- MatLab of mine?
> Leon
I don't know. I do not have MATLAB linprog and I try to
simulate it with my own LP solver obviously without succeed.
What I want to solve is this continue pb
minimize cost * x with constraints
L <= x <= U and a'*x = b
I'm not sure what is the correct syntax of MATLAB linprog to
ignore linear-inequality constraint, I put empty matrix but
MATLAB obviously doesn't like. Hope someone can help.
Otherwise could you try to replace with the following line
with a dummy inequality matrix:
sol = linprog(zeros(size(a')), 1, a', b, cost, L, U);
Bruno
Bruno Luong
>
> sol = linprog(zeros(size(a')), 1, a', b, cost, L, U);
>
My bad, the cost function should be the first argument:
sol = linprog(cost', [], [], a', b, L, U);
OR dummy syntax
sol = linprog(cost', zeros(size(a')), 1, a', b, L, U);
Bruno
Leon galushko
Bruno,
with your last syntax it did!
Now would you know how to get alter the syntax to get
results where (a1,a2, etc) > 0
Leon
Bruno Luong
>
> Now would you know how to get alter the syntax to get
> results where (a1,a2, etc) > 0
> Leon
>
>
Sorry, I don't understand your question. a1, a2 are the
known inputs, what you want to do with them?
Bruno
Leon galushko
Sorry, ich meant c>0 not a.
Bruno Luong
> "Bruno Luong" <b.l...@fogale.findmycountry> wrote in
> message <g6qm7f$5mc$1...@fred.mathworks.com>...
>
> Sorry, ich meant c>0 not a.
>
I see. To ways to do it.
1. Brute force: find all solutions with c>=0, then filter
out by find.
2. Neater way: c>0 means cm1:=c-1 >=0.
So solving
(A) a*c = b such that c > 0
is equivalent to solving
(B) a*cm1 = b-sum(a) such that cm1>=0,
Then take c=cm1+1.
Bruno
Leon galushko
Sorry for pestering you with that but i'm unclear
where would [a*cm1 = b-sum(a) and c=cm1+1] be belonging?
If this refers to MainProg in m-file, thenn to what
syntax sequence have to be that(i thought just you know your
own code anyway better than me bein unexperienced in MathLab).
Thank you
Leon
Bruno Luong
>
> Sorry for pestering you with that but i'm unclear
> where would [a*cm1 = b-sum(a) and c=cm1+1] be belonging?
> If this refers to MainProg in m-file, thenn to what
> syntax sequence have to be that
I let it to you to figure it out Leon.
> (i thought just you know your
> own code anyway better than me bein unexperienced in MathLab).
>
I'm sure you can do it on your own. Don't be hesitate to
play with the code.
Start with simple example, for example:
a = [770 105];
b = 4830;
There is three solutions c>=0, But only two c>0.
Bruno
Leon galushko
O.k I'll do it.
I'm trying to store the solvings in the separated txt file
and alter your syntax:
fid = fopen('C:.../file.txt', 'wt');
disp(c);
fprintf(fid, '\nSolutions c:\n\n');
But it doesn't work!
Leon
Leon galushko
i tryed to see through for myself your function
x = intlin(a, b), but have had difficulties by understanding
it. (i'm not trained mathematican, unfortunately).
I intended to adjust it from x>=0 to x>0, but must say that
my knowledge in linear algebra yet lacking for it.
I would like to use your function for some purposes in my
fee work, but need it really with restrictions there x are
only not null. I would be lucky when you could adjust it for.
Kindly,
Leon