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

发布于 2024-09-25 01:36:31 字数 692 浏览 4 评论 0原文

我有一些关于 {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 和 x^2+y^2<=200 只是示例;实际上,fg 可以是任意次数的任意多项式函数。

编辑:也欢迎 C# 代码。

I have a few inequalities regarding {x,y}, that satisfies the following equations:

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

Note that x and y must be integer.

Graphically it can be represented as follows, the blue region is the region that satisfies the above inequalities:

alt text

The question now is, is there any function in Matlab that finds every admissible pair of {x,y}? If there is an algorithm to do this kind of thing I would be glad to hear about it as well.

Of course, one approach we can always use is brute force approach where we test every possible combination of {x,y} to see whether the inequalities are satisfied. But this is the last resort, because it's time consuming. I'm looking for a clever algorithm that does this, or in the best case, an existing library that I can use straight-away.

The x^2+y^2>=100 and x^2+y^2<=200 are just examples; in reality f and g can be any polynomial functions of any degree.

Edit: C# code are welcomed as well.

如果你对这篇内容有疑问,欢迎到本站社区发帖提问 参与讨论,获取更多帮助,或者扫码二维码加入 Web 技术交流群。

扫码二维码加入Web技术交流群

发布评论

需要 登录 才能够评论, 你可以免费 注册 一个本站的账号。

评论(1

痕至 2024-10-02 01:36:32

对于一般的多项式不等式集,通过枚举搜索以外的任何方法,这当然是不可能的,即使存在有限数量的解。 (也许我应该说这不是微不足道的,因为这是可能的。枚举搜索将起作用,但会受到浮点问题的影响。)请注意,对于高阶不等式,感兴趣的域不需要简单地连接。

编辑:OP 询问如何进行搜索。

问题

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 位机器上运行这个程序,有 8 GB 内存。但即便如此,测试本身仍然非常消耗 CPU。

Elapsed time is 50.182385 seconds.

请注意,浮点考虑有时会导致找到不同数量的点,具体取决于计算的完成方式。

最后,如果您的约束方程更复杂,您可能需要在 y 上的边界表达式中使用根,以帮助确定满足约束的位置。这里的好处是它仍然适用于更复杂的多项式界限。

This is surely not possible to do in general for a general set of polynomial inequalities, by any method other than enumerative search, even if there are a finite number of solutions. (Perhaps I should say not trivial, as it is possible. Enumerative search will work, subject to floating point issues.) Note that the domain of interest need not be simply connected for higher order inequalities.

Edit: The OP has asked about how one might proceed to do a search.

Consider the problem

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

x >= 0, y >= 0

Solve for all integer solutions of this system. Note that integer programming in ANY form will not suffice here, since ALL integer solutions are requested.

Use of meshgrid here would force us to look at points in the domain (0:10000)X(0:10000). So it would force us to sample a set of 1e8 points, testing every point to see if they satisfy the constraints.

A simple loop can potentially be more efficient than that, although it will still require some effort.

% 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

The time required was...

Elapsed time is 0.600419 seconds.

Of the 100020001 combinations that we might have tested for, how many solutions did we find?

size(xy)
ans =
           2     4371264

Admittedly, the exhaustive search is simpler to write.

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

I ran this on a 64 bit machine, with 8 gig of ram. But even so the test itself was a CPU hog.

Elapsed time is 50.182385 seconds.

Note that floating point considerations will sometimes cause a different number of points to be found, depending on how the computations are done.

Finally, if your constraint equations are more complex, you might need to use roots in the expression for the bounds on y, to help identify where the constraints are satisfied. The nice thing here is it still works for more complicated polynomial bounds.

~没有更多了~
我们使用 Cookies 和其他技术来定制您的体验包括您的登录状态等。通过阅读我们的 隐私政策 了解更多相关信息。 单击 接受 或继续使用网站,即表示您同意使用 Cookies 和您的相关数据。
原文