Many design applications can be formulated as optimization constrained by conservation laws. Such optimization can be efficiently solved by the adjoint method, which computes the gradient of the objective to the design variables. Traditionally, the adjoint method has not been able to be implemented in "gray-box" conservation law simulations. In gray-box simulations, the analytical and numerical form of the conservation law is unknown, but the full solution of relevant flow quantities is available. Optimization constrained by gray-box simulations can be challenging for high-dimensional design because the gradient is not efficiently computable.

We consider the case where the flux function is unknown in the gray-box conservation law. The twin model method is presented to estimate the gradient by inferring the flux function from the space-time solution. The method enables the estimation of the gradient by solving the adjoint equation associated with the inferred conservation law. Building upon previous research, a Gaussian process optimization (GPO) framework, the twin-model GPO,  is presented that admits the estimated gradient. The effectiveness of twin-model GPO is compared to a conventional GPO approach where the gradient is unavailable. The performance of conventional GPO approach is found to deteriorate as the optimization dimensionality increases. Adopting the twin-model GPO results in improved objective function evaluation given a limited number of gray-box simulations.