To optimize f(x, y) on a region: find critical points, classify them, evaluate f at critical points and on the boundary, then compare. Global extrema on closed bounded sets are guaranteed to exist.
From single-variable calculus, you know the procedure for finding extrema on a closed interval: set the derivative to zero, find interior critical points, evaluate at the endpoints, and compare all values. Multivariable optimization follows exactly the same logic — it just has more geometry. The "interior" is now the open region inside D, the "endpoints" become the boundary curve (or surface) of D, and "derivative = 0" becomes ∇f = 0.
The gradient ∇f = (∂f/∂x, ∂f/∂y) must vanish at any interior local extremum. This gives a system of two equations in two unknowns — typically nonlinear and potentially with multiple solutions. From second-partials-test-hessian you can classify each critical point by computing the Hessian matrix H and its determinant D = f_xx f_yy − (f_xy)². If D > 0 and f_xx > 0, the point is a local min; if D > 0 and f_xx < 0, a local max; if D < 0, a saddle point. Saddle points are the distinctly multivariable phenomenon — f increases in some directions and decreases in others, so the point is neither a local max nor a local min. They have no single-variable analogue (inflection points with f' = 0 are different).
The critical insight for global optimization is the Extreme Value Theorem: a continuous function on a closed bounded region D always attains a global maximum and minimum, and these must occur either at interior critical points or on the boundary of D. The boundary of a 2D region is typically made up of curves — line segments, circular arcs, or other parameterizable pieces. Optimizing f on each boundary piece reduces to a single-variable problem: parameterize the curve (e.g., x = cos t, y = sin t for the unit circle), substitute into f, and use single-variable calculus to find the extrema of the resulting function of t.
After collecting all candidates — values of f at interior critical points plus values at boundary extrema — you simply compare them. The largest is the global max; the smallest is the global min. There is no shortcut: a negative-definite Hessian at an interior point tells you only that it is a local max, not the global max. The boundary can easily produce higher values. This exhaustive comparison is the entire algorithm, and it always works on closed bounded domains.
This framework is the foundation for Lagrange multipliers, which handle equality-constrained optimization more elegantly. Instead of parameterizing the constraint curve and substituting, Lagrange multipliers directly find the points where ∇f is proportional to ∇g (the constraint gradient). The geometric meaning is the same: at a constrained optimum, f cannot increase in any direction tangent to the constraint, which is precisely the condition ∇f = λ∇g. Lagrange multipliers are the multivariable analogue of the endpoint condition in single-variable calculus.