We present a thin‐sheet approach to modeling deformation at diffuse plate boundaries. We specify plate motion boundary conditions, geometry of the boundary, locations and weakness values of faults, the distribution of sheet weakness, and force potentials that result from surface loads, gravitational potential energy, and flow beneath the sheet. Model optimization is achieved by repeating three steps: (1) a work rate functional is minimized to determine a physical model defined by a priori information; (2) an a posteriori model with identical dimensions and parameters, but different functional form, is constructed, to find a better kinematic fit to observations; and (3) a priori parameters are adjusted, based on comparison of a priori and a posteriori models, and subject to considerations such as constitutive relationships, along‐fault slip distribution profiles, or spatial models of weakness based on other data. Diverse observations of velocity, strain rate, fault slip rate, or fault slip direction are weighted and simultaneously compared with predictions, without bias introduced by arbitrary grouping. We propose an efficient method for estimating optimized model parameters (sheet and fault weakness values and force potentials) that has a wide range of applications: synthesis of observations and boundary conditions into a uniform kinematic model of plate boundary deformation at Earth's surface, physical consideration of long‐term average weakness/stress at the top of the lithosphere and fault weakness/traction, and construction of hypotheses to guide future theoretical or observational studies. We illustrate application of the modeling method to plate boundary deformation in northern South Island, New Zealand.