找到满足不等式约束条件的{x,y}离散对

5

我有几个关于{x,y}的不等式,满足以下方程:

x>=0
y>=0
f(x,y)=x^2+y^2>=100
g(x,y)=x^2+y^2<=200

请注意,xy必须是整数。
图形上可以表示为蓝色区域是满足上述不等式的区域: alt text 现在的问题是,Matlab中是否有任何函数可以找到每个可接受的{x,y}对?如果有一种算法可以做到这种事情,我也很高兴听到它。
当然,我们始终可以使用暴力方法来测试每个可能的{x,y}组合,以查看是否满足不等式。但这是最后的手段,因为它很耗时间。我正在寻找一个聪明的算法来完成这项工作,或者在最好的情况下,我可以直接使用现有的库。 x^2+y^2>=100 and x^2+y^2<=200 只是示例;实际上,fg可以是任何次数的多项式函数。
编辑:欢迎C#代码。
1个回答

4

一般情况下,对于任意一组多项式不等式,除了枚举搜索之外,没有其他方法可以解决这个问题,即使有有限数量的解也是如此。(或许我应该说不是不可能,因为枚举搜索可以解决这个问题,但需要注意浮点数问题。) 需要注意的是,对于高阶不等式,感兴趣的定义域不一定是单连通的。

编辑:原作者询问如何进行搜索。

考虑以下问题:

x^3 + y^3 >= 1e12
x^4 + y^4 <= 1e16

x >= 0, y >= 0

解决此系统的所有整数解。请注意,任何形式的整数规划都无法满足此处请求的所有整数解。

在这里使用网格会强制我们查看域(0:10000)X(0:10000)中的点。因此,它将强制我们对一组1e8个点进行采样,并测试每个点是否满足约束条件。

简单的循环可能比那更有效,尽管仍需要一些努力。

% Note that I will store these in a cell array,
% since I cannot preallocate the results.
tic
xmax = 10000;
xy = cell(1,xmax);
for x = 0:xmax
  % solve for y, given x. This requires us to
  % solve for those values of y such that
  %   y^3 >= 1e12 - x.^3
  %   y^4 <= 1e16 - x.^4
  % These are simple expressions to solve for.
  y = ceil((1e12 - x.^3).^(1/3)):floor((1e16 - x.^4).^0.25);
  n = numel(y);
  if n > 0
    xy{x+1} = [repmat(x,1,n);y];
  end
end
% flatten the cell array
xy = cell2mat(xy);
toc

所需时间为...
Elapsed time is 0.600419 seconds.

我们可能测试了100020001种组合,我们找到了多少解决方案?
size(xy)
ans =
           2     4371264

不可否认,穷举搜索更容易编写。

tic
[x,y] = meshgrid(0:10000);
k = (x.^3 + y.^3 >= 1e12) & (x.^4 + y.^4 <= 1e16);
xy = [x(k),y(k)];
toc

我在一台64位机器上运行了这个测试,拥有8GB的内存。但即便如此,这个测试本身也是一个CPU占用量很高的程序。

Elapsed time is 50.182385 seconds.

请注意,浮点数计算有时会导致找到不同数量的点,具体取决于计算方法。
最后,如果您的约束方程更加复杂,您可能需要在y的边界表达式中使用根来帮助确定约束满足的位置。好处在于它仍然适用于更复杂的多项式边界。

@Ngu:我认为一个问题不可能的一个影响是不能有算法来解决它。我同意@woodchips的观点。如果你对此有疑问,请开始用笔和纸写下前几个满足你的不等式的 (x,y) 对。然后停下来思考一会儿。 - High Performance Mark
对于一般问题,没有简单的解决方案。你可能只能在一个变量上设置一个循环,比如x。将x固定为某个整数值,然后解出满足约束条件的所有y。问题是,要花费一些精力才能证明一旦你到达一个点,在给定x的情况下没有y的解决方案,那么就没有更大的x值可能会产生进一步的解决方案。区间算术可能有助于解决这个问题。 - user85109
@High,如果多项式已知呢?它仍然无法解决吗? - Graviton
仅仅知道多项式系数并不意味着问题变得简单。即使在线性约束的情况下,枚举所有整数解也并非易事。我并不是说它无法解决,而是你会被困在一个循环中。同样地,暴力搜索可以被视为对x的循环,然后对给定x值的所有可能整数解y进行求解。最大的问题在于如何确定何时停止搜索。这对于高阶多项式来说将是困难的。 - user85109
1
@Jouni:不是。按照问题的要求来阅读。原帖中要求找到所有解决方案。线性规划只会得出一个解决方案,而不能找到全部解决方案。 - user85109
显示剩余7条评论

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接