利用OpenGL实现线性规划寻找目标函数极值可视化

简介

最近一段时间学完了简单的线性规划，今天放假闲着于是就花了一点时间把之前用OpenGL实现的函数图像绘制修改成了线性规划的绘图，并且能实现目标函数极值寻找的可视化。

如图这是实现了寻找 z = x^2 + y^2 的最大值的效果

这是实现了寻找 z = 2x + y 的最大值。

当然，非线性的规划也是很容易的，这是寻找 z = x^3 + y^2 的最大值。

说明

Func函数定义你要绘制的二元不等式组，Target函数定义你要寻找极值的目标函数，Color函数一般不用修改，具体我没时间细说了。

代码

#ifdef GL_FRAGMENT_PRECISION_HIGH
precision highp float;
#else
precision mediump float;
#endif

uniform vec2 resolution;
uniform vec3 pointers[10];
uniform int pointerCount;
uniform float time;

float pi = atan(1.0) * 4.0;

//Implicit / f(x) plotter thing.

//Line thickness (in pixels).
#define LINE_SIZE 2.0

//Grid line & axis thickness (in pixels).
#define GRID_LINE_SIZE 1.0
#define GRID_AXIS_SIZE 2.5

//Number of grid lines per unit.
#define GRID_LINES 1.0

const vec2 GRAD_OFFS = vec2(0.005, 0.0);

#define GRAD(f, p) (vec2(f(p) - f(p + GRAD_OFFS.xy), f(p) - f(p + GRAD_OFFS.yx)) / GRAD_OFFS.xx)

//PLOT(Function, Color, Destination, Screen Position)
#define PLOT(f, c, d, p) d = mix(c(d, p), d, smoothstep(0.0, (LINE_SIZE / resolution.y * zoom), abs(f(p) / length(GRAD(f,p)))))

//XY range of the display.
float zoom = 20.0;
vec2 move = vec2(0.0);

float Func(vec2 p)
{
	return float(!(
		(p.x > 0.2) &&
		(p.x + 2.0 * p.y < 10.0) &&
		(2.0 * p.x - p.y < 10.0)
	));
}

float Target(vec2 p)
{
	return p.x * p.x + p.y * p.y;
}

vec3 Color(vec3 d, vec2 p)
{
	return d * 0.8 * (atan(10.0 + -abs(tan(time))
		* Target(p)) / pi + 0.5) + vec3(1, 0, 1);
}

float grid(vec2 p)
{
	vec2 uv = mod(p,1.0 / GRID_LINES);

	float halfScale = 1.0 / GRID_LINES / 2.0;

	float gridRad = (GRID_LINE_SIZE / resolution.y) * zoom;
	float grid = halfScale - max(abs(uv.x - halfScale), abs(uv.y - halfScale));
	grid = smoothstep(0.0, gridRad, grid);

	float axisRad = (GRID_AXIS_SIZE / resolution.y) * zoom;
	float axis = min(abs(p.x), abs(p.y));
	axis = smoothstep(axisRad - 0.05, axisRad, axis);

	return min(grid, axis);
}


void main(void)
{
	vec2 aspect = resolution.xy
		/ resolution.y / 2.0;

	if (pointerCount < 1) {
		move = resolution.y * aspect;
	} else if (pointerCount == 1) {
		move = pointers[0].xy;
	} else if (pointerCount > 1) {
		zoom *= resolution.y / 4.0 /
			distance(pointers[0], pointers[1]);
		for (int n = 0; n < pointerCount; ++n) {
			move += pointers[n].xy;
		}
		move = move / float(pointerCount);
	}

	vec2 uv = ( (gl_FragCoord.xy - move)
		/ resolution.y );
	uv *= zoom;

	vec3 col = vec3(grid(uv) * 0.25 + 0.75);

    PLOT(Func, Color, col, uv);

	gl_FragColor = vec4(vec3(col), 1.0 );
}