Engine map fitting optimization method, device, equipment and storage medium

Abstract

The application relates to an engine MAP fitting optimization method, device, equipment and storage medium. The method comprises the following steps: constructing a target function, performing one-time forward propagation and one-time reverse propagation on the target function through an automatic differentiation framework based on a pre-constructed calculation graph, obtaining a loss function value of the target function and an analytical gradient vector of the target function to all MAP grid node parameters, performing iterative optimization on the MAP grid node parameters by using an analytical gradient-based boundary constraint optimization algorithm based on the loss function value and the analytical gradient vector, and in each iteration, the parameter components reaching the boundary are constrained to the boundary value by using a boundary active set management, parameter updating is performed in the free parameter subspace, and the optimized MAP grid node parameters are obtained. The above method obtains accurate analytical gradients through automatic differentiation of the calculation graph, and combines original projection boundary constraint optimization, solves the gradient distortion and boundary processing problems in high-dimensional MAP calibration, and finally outputs an optimized MAP table meeting engineering requirements.

Description

Technical Field

[0001] This invention relates to the field of calibration and parameter optimization technology for internal combustion engine electronic control units, specifically to an engine MAP fitting optimization method, apparatus, device, and storage medium. Background Technology

[0002] Engine electronic control unit (ECU) calibration is a critical step in the vehicle development process, directly determining the engine's power, fuel economy, and emissions performance. The ECU contains numerous two-dimensional lookup table parameters defined by map plots, such as charging efficiency (MAP), ignition advance angle (MAP), and exhaust gas recirculation (EGR) rate (MAP). These MAP plots typically have a 20x20 or 5x5 grid resolution, with hundreds to thousands of parameters. Calibration engineers need to use numerical optimization methods to fit the grid node values ​​of these MAP plots based on measured data from thousands of steady-state operating points collected from bench tests, minimizing the root mean square error (RMSE) between the predicted values ​​calculated based on MAP interpolation and the measured values. Unlike traditional calibration tasks where each MAP plot is treated as an independent lookup table, modern engine control strategies often combine multiple MAP plots through algebraic expressions to form composite control variables. Therefore, modern calibration tasks present a significant challenge.

[0003] In existing technologies, evolutionary algorithms based on random search can typically find the optimal solution by randomly exploring the parameter space through a population, while online schemes based on iterative learning control can correct the MAP diagram online through closed-loop feedback. However, existing methods generally lack flexible support for user-defined expressions and cannot guarantee a smooth transition of the surrounding area after the locking point is fixed, resulting in insufficient physical plausibility of the fitting results.

[0004] Therefore, how to effectively solve the high-dimensional nonlinear optimization problem of engine MAP fitting is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0005] One objective of this invention is to provide a method, apparatus, device, and storage medium for engine MAP fitting optimization, in order to solve the high-dimensional nonlinear optimization problem in engine MAP fitting.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0007] An engine MAP fitting optimization method includes:

[0008] Construct an objective function, which includes an error term and at least one penalty term for constraining the characteristics of the MAP parameters;

[0009] Based on a pre-built computational graph, the objective function is subjected to one forward propagation and one backward propagation through an automatic differentiation framework to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters.

[0010] Based on the loss function value and the analysis gradient vector, the MAP grid node parameters are iteratively optimized using an analysis gradient-based bounded constraint optimization algorithm. In each iteration, the parameter components that touch the boundary are constrained to the boundary value by the boundary active set management, and the parameter update is performed in the free parameter subspace to obtain the optimized MAP grid node parameters.

[0011] Based on the aforementioned techniques, this method avoids boundary truncation errors and matrix singularity issues associated with finite difference calculus by accurately analyzing gradients. Simultaneously, it utilizes boundary active set management to constrain parameter components touching the boundary to boundary values ​​and performs parameter updates within the free parameter subspace. This maintains the effectiveness of the search direction under tight boundary constraints, achieving stable convergence of high-dimensional MAP parameter optimization under high-dimensional tight boundary constraints.

[0012] Furthermore, the boundary-constrained optimization algorithm based on analytical gradients is a boundary-constrained quasi-Newton optimization algorithm; the parameter update in the free parameter subspace is performed as quasi-Newton optimization in the free parameter subspace.

[0013] Based on the above technical means, a quasi-Newton optimization algorithm with boundary constraints is adopted to perform quasi-Newton optimization in the free parameter subspace. Hessian approximation is constructed using historical gradient information to accelerate convergence. Compared with the first-order gradient method, it significantly reduces the number of iterations at the same accuracy. At the same time, it is combined with boundary active set management to accurately constrain boundary parameters, so that the optimization process can stably satisfy the boundary constraint conditions while maintaining a high convergence speed.

[0014] Furthermore, based on the pre-constructed computational graph, the objective function is subjected to one forward propagation and one backward propagation through an automatic differentiation framework to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters, including:

[0015] Using the MAP grid node parameters as input, forward propagation is performed along the computation graph to obtain the loss function value;

[0016] Perform a backpropagation along the computation graph, and use the chain rule to propagate the gradient of the loss function value back layer by layer to obtain the analysis gradient vector.

[0017] Based on the above technical means, the loss function value is obtained by performing forward propagation on the computation graph, and the analytical gradient vector is obtained by using the chain rule for back propagation. The complete gradient information required for optimization can be obtained simultaneously with one forward and one back propagation, avoiding multiple function evaluations in the finite difference method.

[0018] Furthermore, the penalty term also includes a concentric circle transition penalty term for the locking point;

[0019] The concentric circle transition penalty term for the locking point divides the neighborhood into multiple rings based on the Manhattan distance from each locking point. Neighborhood points in different rings are assigned different attenuation weights, and the attenuation weights decrease from the inside out in each ring, so that the closer the neighboring point is to the locking point, the greater the penalty it receives. This is used to apply smooth transition constraints to the free parameters around the user-selected locking grid point.

[0020] Based on the above technical means, by dividing the neighborhood into multiple rings with Manhattan distance as the center of the locking point and setting the attenuation weight that decreases from the inside to the outside, the free parameters around the locking point are smoothly and asymptotically approached to the locking value, thus avoiding the step-like discontinuity caused by the locking point being fixed alone.

[0021] Furthermore, the method also includes:

[0022] When constructing the parameter boundaries of the objective function, in response to the locked grid points and locked values ​​selected by the user on the MAP graph, the upper and lower bounds of the parameter elements corresponding to the locked grid points are set within the tolerance range of the locked values, so as to apply an immovable hard constraint to the locked grid points.

[0023] Based on the above technical means, by setting the upper and lower bounds of the lock point parameter to within the lock value tolerance range, the mobility of the lock point is directly deprived from the boundary level of the optimizer parameter, making the lock point absolutely immovable during the optimization process.

[0024] Furthermore, the method also includes:

[0025] During the iterative optimization of the MAP grid node parameters, after each parameter update and when the final result is output, the parameter value corresponding to the user-selected locked grid point is forcibly restored to the locked value.

[0026] Based on the above technical means, by forcibly restoring the lock point value after each parameter update and at the final output, lock point drift caused by floating-point precision accumulation or abnormality is eliminated from the software implementation level.

[0027] Furthermore, the penalty term includes a smoothing penalty term;

[0028] The smoothing penalty term is composed of the weighted sum of squared second-order differences in the X direction, the weighted sum of squared second-order differences in the Y direction, and a cross term, and is used to constrain the smoothness of the MAP surface.

[0029] Based on the aforementioned techniques, local curvature is suppressed by weighted summation of second-order differences in the X and Y directions. Smooth surfaces incurves are penalized less, while peaks / concaves are penalized more. Cross terms penalize bowl-shaped / dome-shaped curvature, encouraging flat or linear variations while preserving the true saddle-shaped characteristics of the engine map. This prevents drastic fluctuations in the MAP due to overfitting without sacrificing fitting accuracy.

[0030] Furthermore, the penalty term also includes a gradient constraint penalty term;

[0031] The gradient constraint penalty term applies a unidirectional penalty to the first-order difference that exceeds the preset upper and lower limits of the gradient, which is used to constrain the rate of change of MAP in the coordinate axis direction.

[0032] Based on the above technical means, by applying a unidirectional quadratic penalty to the first-order difference that exceeds the preset upper and lower limits of the gradient, the rate of change of MAP in the X / Y direction is limited to the physically permissible range, thus avoiding violation of the physical laws of the engine.

[0033] Furthermore, based on the loss function value and the analytical gradient vector, an analytical gradient-based boundary-constrained optimization algorithm is used to iteratively optimize the MAP grid node parameters. In each iteration, boundary active set management is used to constrain parameter components touching the boundary to the boundary values, and parameter updates are performed in the free parameter subspace to obtain the optimized MAP grid node parameters, including:

[0034] Step a: In the current iteration step, determine the gradient direction based on the analyzed gradient vector;

[0035] Step b: Based on the gradient direction, in each iteration, the parameter components that touch the boundary are constrained to the boundary value by the boundary active set management, and quasi-Newton optimization is performed in the free parameter subspace to determine the search direction and step size.

[0036] Step c: Update the MAP grid node parameters according to the search direction and step size;

[0037] Repeat steps a-c above until the convergence condition is met to obtain the optimized MAP grid node parameters.

[0038] Based on the above technical means, the boundary-constrained quasi-Newton optimization based on gradient analysis is refined into a specific execution process of "determining gradient direction → boundary active set management to determine search direction and step size → updating parameters → iterating until convergence", so that the optimization algorithm is fully implemented at the claim level.

[0039] Furthermore, the method also includes:

[0040] During the iterative optimization process, the first round starts from the original MAP value, performs a complete iterative optimization, and records the obtained optimal solution as the safety net benchmark.

[0041] In subsequent rounds, starting from the current optimal solution, the MAP grid is downsampled to a low resolution, random noise is applied to the low-resolution grid, and then the original resolution is restored through interpolation upsampling, generating multiple perturbation-based candidate starting points;

[0042] The optimal error is selected from the multiple candidate perturbation starting points by performing forward propagation evaluation, and the actual starting point is then re-executed from the actual starting point for a complete iterative optimization.

[0043] If the results obtained from subsequent iterations are worse than the safety net benchmark and the degree of deterioration exceeds a first preset threshold, then the process reverts to the safety net benchmark and re-perturbs.

[0044] Based on the above technical means, a naturally smooth perturbation MAP is generated by spatially correlated perturbation to avoid non-physical jagged deformation. High-quality perturbation candidates are screened with extremely low computational cost through forward propagation for rapid evaluation. The safety net benchmark prevents cascading failures that lead to "the more you optimize, the worse it gets".

[0045] Furthermore, the method also includes:

[0046] If the error obtained from multiple consecutive optimization rounds does not decrease compared to the previous round, then increase the perturbation rate;

[0047] If the error in this round decreases by more than the second preset threshold compared to the previous round, then the perturbation rate is reduced.

[0048] Based on the above technical means, the perturbation rate is automatically adjusted by monitoring whether the error decreases or improves in multiple consecutive rounds and whether the improvement exceeds the threshold, so as to achieve an automatic balance between "exploration" and "development" without the need for manual preset.

[0049] Furthermore, before constructing the objective function, the method further includes:

[0050] Receive user-inputted custom fitting expression, parse the syntax structure of the expression, extract MAP lookup function calls, data column dependencies, and dependencies between Function nodes, and construct a computational dependency topology in the form of a directed acyclic graph;

[0051] Pre-compute the non-dependent nodes in the computation-dependent topology, and retain the expression strings for the dependent nodes in the computation-dependent topology for subsequent dynamic evaluation;

[0052] By replacing the discontinuous branch selection function in the expression with a continuously differentiable transition function, the computational graph of the traceable gradient propagation path is obtained.

[0053] Based on the above technical means, a differentiable computation graph is constructed by parsing user-defined expressions, so that any composite expression can be automatically differentiated. Combined with dependency analysis and bswitch continuous processing, the calibration tool is generalized, eliminating the need to rewrite and optimize code for each expression.

[0054] Furthermore, the method also includes:

[0055] For each MAP lookup function in the computation graph, the interpolation weight and index of each data point are pre-calculated, and the interpolation weight and index are registered as a fixed buffer for lookup calls during the forward propagation process.

[0056] Based on the above technical means, by pre-calculating the interpolation weights and indices once and registering them as fixed buffers, the interpolation calculation of the MAP lookup function nodes during the forward propagation process is degenerated into a pure lookup multiplication and addition operation. There is no need to repeat the coordinate positioning and interpolation coefficient calculation in each iteration, which significantly reduces the amount of computation per forward propagation. Combined with the finite difference-free characteristic of gradient analysis, the overall optimization efficiency is greatly improved.

[0057] An engine MAP fitting and optimization device, comprising:

[0058] A construction module is used to construct an objective function, which includes an error term and at least one penalty term for constraining the characteristics of the MAP parameters;

[0059] The computation module is used to perform one forward propagation and one backward propagation on the objective function based on a pre-built computation graph through an automatic differentiation framework, to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters;

[0060] The iterative optimization module is used to iteratively optimize the MAP grid node parameters based on the loss function value and the analysis gradient vector using a boundary-constrained optimization algorithm based on the analysis gradient. In each iteration, the parameter components that touch the boundary are constrained to the boundary value by the boundary active set management, and the parameter update is performed in the free parameter subspace to obtain the optimized MAP grid node parameters.

[0061] An electronic device includes: a memory and a processor;

[0062] The memory stores computer-executed instructions;

[0063] The processor executes the computer execution instructions stored in the memory to implement the engine MAP map fitting optimization method described above.

[0064] A computer-readable storage medium storing computer-executable instructions, which, when executed by a processor, are used to implement the engine MAP fitting optimization method as described in any of the preceding claims.

[0065] A computer program product includes a computer program that, when executed by a processor, implements the engine MAP map fitting optimization method described in any of the preceding claims.

[0066] The beneficial effects of this invention are:

[0067] (1) By using an automatic differentiation framework, the loss function value and the accurate analytical gradient vector are obtained simultaneously in one forward propagation, completely avoiding the boundary truncation error of finite difference and the problem of the L-BFGS-B Hessian inverse matrix being destroyed by pseudo-zero gradients. Compared with the optimization stagnation caused by the JTJ matrix singularity of the traditional LM algorithm, the analytical gradient-driven L-BFGS-B scheme is not affected by parameter dependence, fundamentally eliminating the singularity problem. Compared with the catastrophic divergence of finite difference L-BFGS-B, the analytical gradient scheme provides a stable and reliable convergence path.

[0068] (2) L-BFGS-B's native projection gradient mechanism automatically fixes the parameters that are about to touch the boundary to the boundary active set in each iteration step, without the need to manually implement boundary truncation or reflection projection logic. The projection operation is precisely handled by Cauchy point calculation and subspace minimization, which correctly manages the boundary active set, ensuring that the parameters are strictly within the preset range without disrupting the optimizer's search direction.

[0069] (3) Users can customize algebraic fitting expressions of any form through a graphical interface. The system automatically parses and constructs a differentiable computation graph, eliminating the need to write different optimization code for different fitting formulas, which greatly reduces the usage threshold for calibration engineers. The expression framework supports free combination of three parameter types: Map (two-dimensional interpolation parameters), Curve (one-dimensional interpolation parameters), and Scalar (scalar parameters). Combined with recursive node dependency graph construction and dual-path computation strategy, it balances flexibility and computational efficiency.

[0070] (4) Bilinear interpolation pre-computation combined with the finite difference-free gradient analysis avoids n extra forward propagations. The overall acceleration effect reduces the single evaluation time to the millisecond level, and the overall optimization speed is about 40 times faster than the traditional finite difference scheme.

[0071] (5) Spatial correlation perturbation generates a naturally smooth perturbation MAP through a downsampling-upsampling process, avoiding the non-physical jagged deformation produced by traditional independent noise perturbation methods, thus making the perturbation candidates physically reasonable. Rapid multi-candidate screening significantly improves the quality of perturbation candidates with negligible computational cost. Gradient-guided perturbation intensity allocation reduces the number of highly sensitive perturbations and increases the number of low-sensitivity perturbations, further improving perturbation efficiency. Adaptive perturbation rate adjustment enables automatic optimization of the search range. The safety net mechanism effectively prevents the cascading failure of perturbation strategies.

[0072] (6) Under the premise of strictly maintaining the fixed value of the lock point, the free parameters in the area surrounding the lock point are smoothly transitioned to the lock value through the concentric circle transition penalty mechanism with three-circle attenuation weights. The matching forced recovery mechanism eliminates any possible lock point drift at the implementation level. Compared with simple lock point fixation or single global smooth constraint, it achieves the ideal but difficult-to-achieve effect of "smooth transition with fixed point" in engineering.

[0073] (7) The constraint signature inheritance mechanism uses a dual check of "boundary out-of-bounds check + constraint signature comparison" to intelligently determine whether the optimal solution from the previous run still applies to the modified constraints. Any mismatch will trigger a restart from the original imported values, avoiding cross-run parameter pollution. Compared to restarting from the original values ​​every time or directly inheriting without checks, a balance is achieved between safety and efficiency.

[0074] (8) The panel statistical diagnostic system provides a complete analysis chain from macro-statistics to micro-residuals. All panels share a unified vectorized prediction calculation engine, avoiding redundant calculations. Engineers can complete the fitting quality assessment in one stop.

[0075] (9) The end-to-end process from data import to analysis and diagnosis covers the complete business requirements of engine MAP fitting. The entire process is operated through a graphical interface, and users do not need to write code. All parameter values ​​in the user interface are directly transmitted to the engine without any code-level truncation or modification, ensuring that users have complete control over the optimization process. Attached Figure Description

[0076] Figure 1 This is the application system architecture diagram for this solution;

[0077] Figure 2 Example 1 of the engine MAP diagram fitting and optimization method provided for this solution;

[0078] Figure 3 Example 2 of the engine MAP map fitting and optimization method provided for this solution;

[0079] Figure 4 Example 3 of the engine MAP map fitting and optimization method provided for this solution;

[0080] Figure 5 Example 4 of the engine MAP map fitting and optimization method provided for this solution;

[0081] Figure 6 Example 5 of the engine MAP plot fitting and optimization method provided for this solution;

[0082] Figure 7 A schematic diagram of the engine MAP fitting and optimization device provided for this solution;

[0083] Figure 8 A schematic diagram of the electronic device structure provided for this solution. Detailed Implementation

[0084] The embodiments of the present invention will be described below with reference to the accompanying drawings and preferred embodiments. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific embodiments, and various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be understood that the preferred embodiments are only for illustrating the present invention and not for limiting the scope of protection of the present invention.

[0085] It should be noted that the illustrations provided in the following embodiments are only schematic representations of the basic concept of the present invention. Therefore, the drawings only show the components related to the present invention and are not drawn according to the actual number, shape and size of the components in the actual implementation. In the actual implementation, the form, quantity and proportion of each component can be arbitrarily changed, and the layout of the components may also be more complex.

[0086] Figure 1 This is the application system architecture diagram for this solution, such as... Figure 1 As shown in the diagram, the system has a three-layer structure: a graphical user interface layer (based on an open-source graphical user interface framework) that provides the interface for interaction between the user and the system. This layer includes a data import page, a model configuration page, a Function node editing page, a parameter configuration page, an optimization execution page, and an analysis and evaluation page, corresponding to each stage of the complete workflow from data import to analysis and evaluation. Users can complete the input of custom fitting expressions, configuration of MAP parameters, initiation of optimization tasks, and viewing and analysis of results through the graphical interface without writing any code.

[0087] The data model layer carries the data structures and configuration information upon which the entire optimization task depends. This layer includes an engine electronic control unit (ECU) standard description file parser, used to parse the ECU's standard description file to extract MAP definition information; a data mapping module, used to establish the correspondence between measured data and MAP parameters; a MAP parameter structure, used to store MAP grid node values ​​and coordinate axis information; a function node module, used to represent function nodes after parsing user-defined expressions and their dependencies; and a dataset module, used to manage the measured data collected from bench tests.

[0088] The optimization engine layer, implemented based on an automatic differentiation framework and a numerical optimization library, is the computational core of the entire system. This layer includes: a differentiable objective function encapsulation module, responsible for encapsulating the entire optimization problem into an automatically differentiable computational module; a pre-computation bilinear interpolation module, used to pre-compute the interpolation weights and indices of all data points and register them as fixed buffers before optimization begins; a composite loss function module, used to calculate the multi-objective loss function value composed of the sum of squared errors, a smoothing penalty term, a lock-point concentric circle transition penalty term, and a gradient constraint penalty term; and a quasi-Newton method solver interface with boundary constraints, used to receive the loss function value and analyze the gradient vector, performing iterative optimization with projected gradient boundary constraints.

[0089] Figure 2 Example 1 of the engine MAP plot fitting optimization method provided for this solution, as follows: Figure 2 As shown, the method includes:

[0090] S201: Construct the objective function.

[0091] This step constructs a differentiable multi-objective loss function as the driving objective for subsequent optimization iterations. This objective function measures the fitting deviation between the MAP predicted values ​​and the measured values ​​on the test bench through an error term, and constrains the physical rationality of the MAP parameters through a penalty term. This ensures that the optimization process not only pursues fitting accuracy but also guarantees that the MAP surface meets engineering requirements such as smoothness, physical boundaries, and user-locked points. The objective function includes an error term and at least one penalty term used to constrain the characteristics of the MAP parameters.

[0092] Specifically, the mathematical representation of the objective function can be expressed as:

[0093] L = SSE + smoothing penalty + lockpoint concentric circle transition penalty + gradient constraint penalty

[0094] Where SSE is the sum of squared errors between the predicted and measured values.

[0095] Optionally, the penalty term may include any one or more combinations of the lockpoint concentric circle transition penalty term, smoothing penalty term, and gradient constraint penalty term.

[0096] If only fitting accuracy is required, only the SSE error term can be used; if overfitting needs to be prevented, a smoothing penalty term can be added on top of the SSE; if the rate of change of MAP needs to be limited, a gradient constraint penalty term can be added; if a specific grid point needs to be locked and the surrounding area needs to be smoothly transitioned, a concentric circle transition penalty term for the locked point can be added; for complex engineering requirements, all four components can be used simultaneously.

[0097] In this process, using the MAP grid node values ​​as input, bilinear interpolation is performed on the query coordinates of each data point to calculate the predicted value corresponding to that point. The difference between the predicted and measured values ​​is squared, and the squared errors of all data points are summed to obtain the sum of squared errors (SSE). SSE is the dominant term of the objective function, and minimizing it ensures that the MAP surface passes through the measured data points as accurately as possible.

[0098] The smoothing penalty term consists of the weighted sum of squared second-order differences in the X direction, the weighted sum of squared second-order differences in the Y direction, and cross terms, and is used to constrain the smoothness of the MAP surface.

[0099] The second-order gradient squared term is used to ensure that the curvature of the surface changes smoothly. It is achieved by calculating the second-order difference in the X and Y directions and summing the squares. The cross term is used to ensure the coordination between gradients in each direction and to avoid non-physical distortion of the surface.

[0100] Each item is multiplied by its corresponding weight coefficient (the smoothness weight set by the user for each MAP in the interface) and then summed to obtain the smoothing penalty term.

[0101] The gradient constraint penalty term is activated when the user sets finite upper and lower limits for the gradient. This term uses the parameter smoothness weight as a baseline, multiplied by a preset amplification factor (approximately 100 times) as the total weight coefficient, and applies a unidirectional quadratic penalty to the first-order difference that exceeds the user-defined upper and lower limits for the gradient in the X and Y directions.

[0102] Specifically, the system calculates the first-order difference value of the MAP for each pair of adjacent grid points in the X and Y directions, and determines whether the difference value exceeds the user-preset upper or lower limit of the gradient. For the portion exceeding the upper limit, the excess is calculated and a secondary penalty is applied; for the portion below the lower limit, the excess is calculated and a secondary penalty is applied; difference values ​​that do not exceed the limit are not penalized. All excess penalty values ​​are summed and multiplied by the total weight coefficient to obtain the gradient constraint penalty term.

[0103] The purpose of setting the above amplification factor is to ensure that the gradient constraint penalty term has a perceptible proportion in the total loss, so that it is not ignored by the optimizer during the optimization process, thereby effectively exerting the constraint effect.

[0104] The concentric circle transition penalty term for the lock point divides the neighborhood into multiple rings based on the Manhattan distance from each lock point. Neighborhood points in different rings are assigned different attenuation weights, and the attenuation weights decrease from the inside out in each ring. This results in a greater penalty for neighboring points that are closer to the lock point. This is used to apply smooth transition constraints to the free parameters around the user-selected lock grid point.

[0105] This term uses the parameter smoothness weight multiplied by a fixed coefficient as the total weight, and then multiplied by the sum of the squared deviations of the attenuation weights of each cycle.

[0106] For example, with each lock point as the center, three concentric rings of neighborhood are formed outwards according to Manhattan distance: the first ring consists of four neighboring points at a distance of 1 (4 in total), the second ring consists of eight neighboring points at a distance of 2 (8 in total), and the third ring consists of twelve neighboring points at a distance of 3 (12 in total). The decay weight configuration for each ring of neighborhood points is as follows: weight w1=1.0 for the first ring, weight w2=0.5 for the second ring, and weight w3=0.25 for the third ring.

[0107] The squared deviation between the parameter values ​​of each concentric circle's neighboring points and the locked point value is calculated, multiplied by the corresponding concentric circle's attenuation weight, and then summed. Finally, it is multiplied by the total weight coefficient to obtain the concentric circle transition penalty term for the locked point. This penalty term makes the constraints stronger for neighboring points closer to the locked point and weaker for neighboring points farther away, thereby achieving a smooth asymptotic transition around the fixed point.

[0108] In one implementation, when constructing the parameter boundaries of the objective function, in response to the user's selection of locked grid points and locked values ​​on the MAP graph, the upper and lower bounds of the parameter elements corresponding to the locked grid points are set to the tolerance range of the locked values, so as to impose an immovable hard constraint on the locked grid points.

[0109] During the iterative optimization of MAP grid node parameters, after each parameter update and when the final result is output, the parameter values ​​corresponding to the user-selected locked grid points are forcibly restored to the locked values.

[0110] When constructing the parameter boundaries of the objective function, in response to the locked grid points and their locked values ​​selected by the user on the MAP graph, the upper and lower bounds of the parameter elements corresponding to the locked grid points are set to the locked values ​​plus or minus the minimum tolerance, so that the parameter elements are immovable at the variable boundary level of the optimizer.

[0111] During the optimization iteration process, a forced recovery operation is performed after each parameter update to restore the parameter values ​​corresponding to the locked grid points to their locked values. This operation is performed after each iteration step parameter update, on all smoothing post-processing paths, and when the final result is written back, forming a triple guarantee.

[0112] The aforementioned hard constraints and forced recovery mechanisms work in conjunction with the concentric circle transition penalty for the lock point: hard constraints ensure that the lock point remains absolutely stationary, the concentric circle transition penalty ensures that the area around the lock point is smoothly asymptotically approached, and forced recovery eliminates any possible lock point drift from the software implementation level. Together, these three constitute a complete lock point constraint coding scheme.

[0113] S202: Based on a pre-built computational graph, the objective function is propagated once forward and once backward through an automatic differentiation framework to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters.

[0114] This step is based on a pre-built computation graph. It performs forward propagation on the objective function using an automatic differentiation framework to calculate the loss function value, and then uses the chain rule to calculate the gradient value of the loss function with respect to the parameters of each MAP grid node through back propagation, providing a precise search direction for subsequent gradient-driven optimization.

[0115] Specifically, using the MAP grid node parameters as input, forward propagation is performed along the computation graph to obtain the loss function value. Then, backward propagation is performed along the computation graph, and the gradient of the loss function value is passed backward layer by layer using the chain rule to obtain the analysis gradient vector.

[0116] S203: Based on the loss function value and the analytical gradient vector, the MAP grid node parameters are iteratively optimized using an analytical gradient-based optimization algorithm with boundary constraints. In each iteration, the parameter components that touch the boundary are constrained to the boundary value by the boundary active set management. The parameter update is performed in the free parameter subspace to obtain the optimized MAP grid node parameters.

[0117] This step uses the precise analytical gradient vector and loss function values ​​obtained in S202 as driving information to iteratively update the MAP grid node parameters through an optimization algorithm with boundary constraints. In each iteration, the boundary active set management mechanism is used to identify parameter components that have reached the boundary and fix them to the boundary values. Parameter updates are performed only in the free parameter subspace, enabling the optimization to converge efficiently under strict boundary constraints. Finally, the optimal MAP grid node parameters that satisfy all constraints are output.

[0118] The core of this step is a general boundary constraint processing mechanism, whose workflow can be summarized into four stages: identification, fixation, dimensionality reduction, and update.

[0119] In the current iteration step, boundary detection is first performed on each component of the parameter vector to determine whether it touches a preset lower or upper bound. Parameter components that are determined to have touched the boundary are marked as "active", forming the current active boundary set; parameter components that have not touched the boundary are marked as "free variables".

[0120] After labeling, all active parameter components are fixed to their current boundary values, preventing them from participating in the update calculations of this iteration. Simultaneously, these active parameters are removed from the optimization subspace, leaving only the free variables to form the free parameter subspace. Through this dimensionality reduction, optimization calculations are performed only within the free parameter subspace, while the boundary parameters remain unchanged.

[0121] Within the free parameter subspace, the search direction is determined using the components of the corresponding free parameters in the current analysis gradient vector. After determining the step size through line search, updates are performed only on the free parameter components. Parameters in the boundary active set remain unchanged throughout the entire process.

[0122] It's important to note that the aforementioned "boundary active set management + free parameter subspace optimization" mechanism is independent of the specific optimization algorithm used. Whether it's gradient descent, conjugate gradient, quasi-Newton methods, or trust region methods, this mechanism can handle boundary constraints. The only difference lies in the specific strategies for determining the search direction and step size within the free parameter subspace. However, all strategies are based on analyzing gradient vectors to drive iterations and using active set management to precisely constrain parameters touching the boundary values. When the current optimal solution lies on the boundary, the boundary parameters remain unchanged, and the free parameters continue to optimize. If a free parameter touches the boundary in a subsequent iteration, it will be identified by the active set and included in the boundary constraints in the next iteration. This dynamically updated mechanism ensures that boundary management is always synchronized with the current iteration state.

[0123] Specifically, the current parameter vector and its corresponding analytical gradient vector are used as inputs for each iteration. Each component of the current parameter vector has preset physical boundary constraints, namely, a lower bound and an upper bound.

[0124] In the current iteration step, boundary checks are performed on each component of the parameter vector. If the value of a component is less than or equal to its lower bound (within a preset minimum tolerance range), the component is marked as "lower bound active" and fixed at the lower bound value; if the value of a component is greater than or equal to its upper bound (within a preset minimum tolerance range), the component is marked as "upper bound active" and fixed at the upper bound value; if a component is between the upper and lower bounds and far from the boundary, it is marked as a "free variable".

[0125] All parameter components marked as active constitute the current boundary active set.

[0126] Once the boundary active set is determined, these active parameter components are removed from the optimization subspace. Specifically, in this iteration, the active parameter components do not participate in parameter update calculations, and their values ​​remain unchanged from the boundary values. Only the parameter components marked as "free variables" and their corresponding gradient components are retained, forming the free parameter subspace.

[0127] Within the free parameter subspace, based on the components of the corresponding free parameters in the current analysis gradient vector, any gradient-driven optimization algorithm is used to determine the search direction and update the parameters. The specific strategies of different algorithms within the free subspace are as follows:

[0128] Gradient descent directly uses the negative gradient direction as the search direction and performs a line search along this direction to determine the step size. This method utilizes first-order gradient information, is computationally simple, and has a small computational cost per step, making it suitable for scenarios requiring rapid acquisition of feasible solutions.

[0129] The conjugate gradient method constructs a conjugate direction by combining historical search directions with the current gradient, and then performs a line search along this conjugate direction to determine the step size. This method accelerates convergence compared to gradient descent, is suitable for medium-sized optimization problems, and can achieve an exact solution to a quadratic function within a finite number of steps.

[0130] The quasi-Newton method constructs an approximation of the Hessian matrix using historical gradient information, determines the search direction by solving the approximate Newton equation, and then performs a line search along this direction to determine the step size. This method utilizes second-order curvature information to accelerate convergence and is one of the mainstream methods for boundary-constrained optimization problems. It has a fast convergence speed and is suitable for high-dimensional optimization problems.

[0131] The projective conjugate gradient method adds a projection operation to the conjugate gradient method. After updating the parameters along the conjugate direction in the free parameter subspace, out-of-bounds parameters are projected. This method combines the fast convergence of the conjugate gradient method with the boundary handling capability of the projection operation, making it suitable for medium-sized optimization problems with many boundary constraints.

[0132] The trust region method constructs an approximate model of the objective function in each iteration, solves a subproblem with boundary constraints within the trust region radius, and determines whether to accept an iteration step and adjusts the trust region radius by comparing the actual descent with the predicted descent. This method has good theoretical guarantees of global convergence, and the solution of the subproblems typically relies on the conjugate gradient method or other iterative methods.

[0133] Sequential quadratic programming determines the search direction in each iteration by solving a quadratic programming subproblem. This method utilizes the gradient information of the current iteration point and the Hessian approximation to construct a quadratic programming subproblem. The constraints of this subproblem include the linearized boundary constraints of the original problem. Solving this subproblem yields the search direction, and then the step size is determined through line search. Sequential quadratic programming has a natural advantage in handling boundary constraints, as each iteration's subproblem contains precise constraints on the active boundary variables.

[0134] The algorithms described above are merely a few examples. The method provided in this embodiment is applicable to, but not limited to, the algorithms mentioned above. All of these algorithms share the common characteristics of "iteration driven by analytical gradients, precise handling of boundary constraints through active set management, and parameter updates performed only within the free parameter subspace." The only difference lies in the specific strategies for determining the search direction and step size within the free parameter subspace; however, all strategies are based on analytical gradients and all utilize the same set of boundary active set management mechanisms. Therefore, a suitable implementation method can be selected from the above algorithms or other applicable algorithms based on the scale, nonlinearity, and computational resources of the specific problem.

[0135] After each iteration, the convergence criteria are checked. Common convergence criteria include: whether the maximum component of the projected gradient is less than a preset tolerance, whether the change in the loss function value is less than a preset threshold, whether the change in the parameter vector is less than a preset threshold, and whether the number of iterations has reached a preset upper limit.

[0136] If any convergence condition is met, the iteration terminates, and the current parameter vector is output as the optimal solution. If the convergence condition is not met, the next iteration step is initiated: the loss function value under the current parameters is recalculated and the gradient vector is analyzed; boundary detection is re-executed to update the active set; and the parameters are updated again in the updated free parameter subspace. This process is repeated until convergence is achieved.

[0137] Optionally, the boundary-constrained optimization algorithm based on analytical gradients can be a boundary-constrained quasi-Newton optimization algorithm. Accordingly, performing parameter updates in the free parameter subspace is equivalent to performing quasi-Newton optimization in the free parameter subspace.

[0138] Specifically, when the algorithm is a quasi-Newton optimization algorithm with boundary constraints, the specific optimization process includes: in the current iteration step, determining the gradient direction based on the analyzed gradient vector, using boundary active set management to constrain the parameter components that touch the boundary to the boundary value, performing quasi-Newton optimization in the free parameter subspace, determining the search direction and step size, updating the MAP grid node parameters according to the search direction and step size, repeating the above steps until the convergence condition is met, and obtaining the optimized MAP grid node parameters.

[0139] It's important to note that in optimization problems with boundary constraints, variables are limited to their respective upper and lower bounds. During the iteration process, some variables will touch the boundary (reaching the upper or lower bound). If these variables "on the boundary" continue to participate in quasi-Newton updates, problems may arise such as exceeding the boundary after the update, requiring additional handling, reduced line search efficiency, distortion of the search direction near the boundary, and variables oscillating back and forth on the boundary, leading to slower convergence. Therefore, a mechanism is needed to manage these variables on the boundary, i.e., boundary active set management.

[0140] In constraint optimization, the active set refers to the set of variables that make a certain inequality constraint equal at the current iteration point. In this embodiment, the active set refers to the set of MAP grid node parameters that have reached the upper or lower bound in the current iteration. These parameters are fixed at the boundary value in this iteration and do not participate in the quasi-Newton update.

[0141] Boundary active set management identifies which MAP grid node parameters have reached their upper or lower bounds in each iteration, removes these parameters from the current optimization subspace and fixes them on the boundary values, so that quasi-Newton updates are performed only on the remaining free parameters, thereby achieving efficient optimization with boundary constraints.

[0142] The engine MAP graph fitting optimization method provided in this application constructs an objective function. Based on a pre-constructed computational graph, it performs one forward propagation and one backward propagation on the objective function through an automatic differentiation framework to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters. Based on the loss function value and the analytical gradient vector, iterative optimization of the MAP grid node parameters is performed using an analytical gradient-based boundary constraint optimization algorithm. In each iteration, boundary active set management is used to constrain the parameter components that touch the boundary to the boundary value. Parameter updates are performed in the free parameter subspace to obtain the optimized MAP grid node parameters. The above method constructs an objective function that includes physical constraint penalties, quantifying engineering experience (such as MAP surface smoothness, rate of change limits, and precise fixation of locking points) into mathematical constraints and integrating them into the optimization objective. Based on this, automatic differentiation is performed based on a pre-constructed differentiable computation graph. The loss function value is calculated with one forward propagation and the gradient vector of the loss function with respect to all parameters to be optimized is accurately obtained with one backpropagation. This avoids the truncation error caused by the difficulty in step size selection in the traditional finite difference method and the computational overhead of multiple function evaluations. At the same time, it fundamentally eliminates the optimization stagnation problem caused by the singularity of the JTJ matrix in the LM algorithm. Furthermore, iterative optimization is performed using a boundary-constrained optimization algorithm based on analytical gradients. In each iteration, the parameter components that have reached the upper and lower bounds are dynamically identified and fixed to the boundary values ​​through boundary active set management. At the same time, gradient-driven updates are performed on the remaining parameters in the free parameter subspace. This boundary handling mechanism ensures that the optimization search direction remains effective within the feasible region, avoiding the problems caused by the penalty function method destroying the objective function shape or simple truncation leading to a distortion of the search direction. Ultimately, it achieves stable and efficient convergence of high-dimensional MAP parameters under the premise of strictly satisfying all physical boundary constraints, and outputs an optimized MAP table that meets engineering requirements.

[0143] Figure 3 Example 2 of the engine MAP fitting optimization method provided for this solution, such as Figure 3 As shown, based on the above embodiment, step S202 specifically includes:

[0144] S301: Using the MAP grid node parameters as input, perform forward propagation along the computation graph to obtain the loss function value.

[0145] This step inputs the current values ​​of the MAP grid node parameters into the computation graph, calculates each node from input to output along the topological order of the computation graph, and finally outputs the loss function value of the objective function, which serves as an evaluation index of the fitting quality of the current parameters and the starting point for backpropagation.

[0146] Specifically, the computation graph was constructed during the preparation phase before S301 was executed. This computation graph, in the form of a directed acyclic graph, records the dependencies between nodes in the user-defined fitting expression, including MAP lookup function nodes, algebraic operation nodes, function nodes, and data column reference nodes. The leaf nodes in the computation graph are MAP grid node parameters, and the root node is the loss function node.

[0147] The values ​​of all current MAP grid node parameters (i.e., the variables to be optimized) are used as inputs to the computation graph, and the assignment is completed at the input layer.

[0148] The computation logic of each node is executed sequentially according to the topological order of the computation graph (from the leaf node to the root node). For each node, after all its dependent predecessor nodes have completed their computations, the node reads the output value of the predecessor node as input, executes its own computation logic, and outputs the computation result to its successor node.

[0149] For a MAP lookup function node, the operation logic is as follows: This node reads the coordinates of the data point, and uses the pre-calculated interpolation weights registered as a fixed buffer. ) and one-dimensional flattened grid index ( The system retrieves the corresponding four corner point parameter values ​​through four index lookups, multiplies them by the corresponding interpolation weights, and then sums them (four multiplications and three additions) to obtain the interpolation result for that data point, which is then output to the successor node. Throughout the forward propagation process, the system only performs table lookups and multiplication-addition operations, without performing any searches or dynamic interpolation coefficient calculations.

[0150] For an algebraic operation node, the operation logic is as follows: the node reads the output value of its predecessor node, performs the corresponding algebraic operation (such as addition, subtraction, multiplication, division, exponentiation, etc.), and outputs the calculation result to the successor node.

[0151] For a Function node, the computation logic is as follows: the node reads the output value of its predecessor node, executes the user-defined function operation, and outputs the calculation result to the successor node. If the Function node is determined to be a non-dependent node in dependency analysis, its calculation result has been pre-calculated in the preparation phase, and the pre-calculated result is directly read in the forward propagation; if it is a dependent node, it is dynamically evaluated in the forward propagation.

[0152] For data column reference nodes, the operation logic is as follows: the node directly reads the corresponding test bench measured data value and outputs it as a constant to the successor node.

[0153] The output value of each node is passed down layer by layer along the edges of the computation graph until all nodes have completed the computation.

[0154] The forward propagation eventually reaches the loss function node. This node reads the predicted values ​​and measured data passed from its predecessor node, calculates the sum of squared errors between them, and simultaneously reads the parameter values ​​of the current MAP grid node, calculating the smoothing penalty, the concentric circle transition penalty for the lock point (if the user has selected a lock point), and the gradient constraint penalty (if the user has set upper and lower gradient limits). Each component is then weighted and summed according to preset weight coefficients, outputting the final loss function value f.

[0155] Once the loss function node completes its calculation and outputs the loss function value f, the forward propagation process ends. The current loss function value f is recorded as the starting point for the S302 backpropagation.

[0156] S302: Perform a backpropagation along the computation graph, using the chain rule to propagate the gradient of the loss function value back layer by layer to obtain the analysis gradient vector.

[0157] This step follows the reverse topological order of the computation graph, using the chain rule to propagate the gradient of the loss function value layer by layer in reverse, and finally calculates the analytical gradient vector of the loss function value with respect to all MAP grid node parameters, which serves as the search direction for gradient-driven optimization.

[0158] Specifically, starting with the loss function value f output by S301, the backpropagation interface of the automatic differentiation framework is called. The automatic differentiation framework initializes the gradient of the loss function value f with respect to itself to 1, and propagates the gradient layer by layer backward along the reverse direction of the computation graph, starting from the loss function node.

[0159] During backpropagation, each node is traversed in reverse order of the computation graph topology (from the root node to the leaf nodes). For the currently traversed node, its inputs include: the upstream gradient propagated from its downstream nodes (i.e., the gradient of the loss function value with respect to the node's output value) and the input value recorded by the node during forward propagation.

[0160] Based on the derivative rule of this node (i.e., the Jacobian matrix of the node's operation with respect to the input), the upstream gradient is multiplied by its local derivative with respect to the input to calculate the gradient of the loss function value with respect to each input variable of this node.

[0161] Taking algebraic operation nodes as an example: if a node performs an addition operation... The gradient from downstream Multiply by respectively and ,get and The propagation is performed in reverse order to its predecessor nodes x and y, respectively; if the node performs a multiplication operation... ,but .

[0162] Taking the MAP lookup table function node as an example: using the interpolation weights recorded in the forward propagation, the gradient from the downstream is multiplied by the interpolation weight coefficients of each corner point and distributed to the corresponding four corner point parameters, thus completing the reverse distribution of the gradient.

[0163] The gradients calculated by each node are propagated back along the edges of the computation graph to their predecessor nodes. Each node receives the sum of gradients returned by all its successor nodes (if a node has multiple successors, the sum of the gradients is used as the final upstream gradient of that node) and continues to propagate upstream.

[0164] The gradient is propagated layer by layer through each layer of the composite function, from the loss function node to the penalty term calculation node, then to the algebra operation node, then to the MAP lookup table function node, until it reaches the leaf node of the computation graph. Finally, the gradient is propagated backward to the leaf node corresponding to the MAP grid node parameter.

[0165] After backpropagation has traversed all nodes in the computational graph, the automatic differentiation framework records the gradient of the loss function with respect to the parameters of each MAP grid node. These gradients together form a complete vector g, the dimension of which is equal to the total number of parameters in all MAP grid nodes (e.g., the sum of the number of grid points in multiple MAP tables, such as 825 dimensions). This vector is the analytical gradient vector of the loss function with respect to all optimization variables.

[0166] The loss function value f (scalar) and the analysis gradient vector g (n-dimensional vector, where n is the total number of parameters) together form the (f, g) tuple, which is directly passed to the optimization algorithm with projected gradient boundary constraints as the complete computation result of a single objective function call. No additional computation is required to start the next iteration of optimization.

[0167] The engine MAP graph fitting optimization method provided in this application uses MAP grid node parameters as input, performs forward propagation along the computation graph to obtain the loss function value, performs one backpropagation along the computation graph, and uses the chain rule to propagate the gradient of the loss function value back layer by layer to obtain the analysis gradient vector. This method performs forward and backpropagation along the computation graph, obtaining the loss function value and the accurate analysis gradient vector in a single function call, providing efficient gradient information without truncation error for optimization.

[0168] Figure 4 Example 3 of the engine MAP fitting optimization method provided for this solution, such as Figure 4 As shown, based on the above embodiments, when the boundary-constrained optimization algorithm based on the analysis gradient is a boundary-constrained quasi-Newton optimization algorithm, step S203 specifically includes:

[0169] Step a: In the current iteration step, determine the gradient direction based on the analyzed gradient vector;

[0170] Step b: Based on the gradient direction, in each iteration, the parameter components that touch the boundary are constrained to the boundary values ​​by the boundary active set management, and quasi-Newton optimization is performed in the free parameter subspace to determine the search direction and step size.

[0171] Step c: Update the MAP grid node parameters according to the search direction and step size;

[0172] Repeat steps a-c above until the convergence condition is met to obtain the optimized MAP grid node parameters.

[0173] Using the precise analytical gradient vector and loss function values ​​obtained from S202, the MAP grid node parameters are iteratively optimized using the L-BFGS-B quasi-Newton optimization algorithm with boundary constraints. By using boundary active set management, the parameter components that touch the boundary are constrained to the boundary values. Quasi-Newton optimization is performed in the free parameter subspace to ensure that the parameters are always within the preset physical boundary range. Finally, the optimal MAP grid node parameters that meet the fitting accuracy and physical constraints are obtained.

[0174] Determine the analysis gradient vector Then, the optimizer performs boundary projection constraint processing, based on the user-preset boundary range. Cauchy point calculations are performed, determining the set of dimensions that touch the boundary based on the temporal order in which each dimension reaches the boundary. This boundary extent is a physical lower limit set by the user for each MAP mesh node during the parameter configuration phase. and physical limit constitute.

[0175] For the L-BFSG-B native projection mechanism, the optimizer performs the following subspace minimization operation:

[0176] Fixed the set of dimension variables that touch the boundary to the boundary value;

[0177] In the subspace formed by the remaining free variables, the search direction (quasi-Newtonian direction) is used. Continue to optimize;

[0178] The optimizer uses a two-loop recursion algorithm (L-BFGS-B) to construct the search direction. Specifically, the optimizer maintains the parameter changes from the most recent m iterations. and gradient change ,in The two-loop recursive algorithm utilizes these and The search direction can be directly calculated from a vector through a series of vector dot products and scaling operations. Furthermore, it eliminates the need to explicitly store the complete n×n Hessian matrix or its inverse (the Hessian matrix has an n×n dimension, where n is the total number of parameters; for MAP calibration problems with hundreds to thousands of parameters, storing and inverting the complete Hessian matrix is ​​computationally infeasible). This characteristic makes the L-BFGS-B algorithm particularly suitable for high-dimensional optimization problems.

[0179] After obtaining candidate points, perform a line search along the search direction to determine the appropriate step size. The goal of line search is to find a step size that allows the objective function value to decrease sufficiently while satisfying the Wolfe condition, i.e., satisfying the sufficient descent condition (Armijo condition) and the curvature condition, ensuring that the step size is neither too large nor too small, and avoiding oscillations while maintaining convergence speed.

[0180] After determining the step size, the optimizer updates the parameters according to the following formula: (Free subspace set) or (Boundary active set).

[0181] This projection operation is the native box bounds constraint handling mechanism of the L-BFGS-B algorithm. The L-BFGS-B algorithm identifies the currently active set of boundary constraints by calculating generalized Cauchy points. During the subspace minimization phase, the activated variables are projected back onto the boundaries, while simultaneously preserving the search direction of free variables (variables not touching the boundaries). This method ensures that the parameters are strictly within the feasible region without disrupting the optimizer's search direction, avoiding the search direction failure problem near the boundaries caused by penalty function methods or simple truncation methods.

[0182] Repeat steps a to c. In each iteration, recalculate the loss function value and analyze the gradient vector based on the current parameter position, construct a new gradient direction, calculate the boundary active set to determine the search direction and step size, update the parameters, and form a complete iterative closed loop.

[0183] After the projection operation is completed in each iteration step, if the user has selected the locked grid points and their locked values ​​on the MAP map, the system performs an additional forced recovery operation after the L-BFGS-B projection, directly assigning the parameter values ​​corresponding to the locked grid points to the locked values. This operation is independent of the optimizer's projection gradient mechanism and serves as a supplementary guarantee at the software implementation level. Together with the projection operation in step c, it constitutes a two-layer boundary control (the native box bounds projection of L-BFGS-B is used to handle the general boundaries of all parameters, while forced recovery is used to specifically handle the precise fixing of the locked points).

[0184] After each iteration, the optimizer checks whether the following convergence conditions are met:

[0185] Whether the gradient projection norm is less than the preset tolerance (pgtol), that is, whether the maximum component of the projection gradient is lower than the threshold;

[0186] Whether the change in the loss function value is less than a preset threshold (ftol);

[0187] Has the number of iterations reached the preset upper limit (maxiter)?

[0188] When any of the above conditions are met, the optimizer terminates the iteration, compares the RMSE of the best-seen parameters with the RMSE of the last iteration point. If the best-seen result is better, the parameters corresponding to the best-seen parameters are used to overwrite the final result. Otherwise, the current optimal parameters are output. If the optimizer throws an exception, an exception rollback mechanism is triggered, and the historical best parameters (best-seen) tracked internally by the optimizer are output as the final result. That is, throughout the entire iteration process, the optimizer always maintains a copy of the parameters that have been visited and minimize the loss function value. When an exception occurs, it automatically rolls back to this copy to avoid the deterioration of the output result due to the final step bounce or numerical collapse.

[0189] The entire optimization process employs a two-stage architecture: In the initial optimization phase (multi-starting point exploration phase), a smaller L-BFGS-B memory vector (maxcor=15) is used, making the algorithm more inclined towards global exploration and avoiding premature convergence to local optima; in the later optimization phase (deep fine-tuning phase), the system uses a larger L-BFGS-B memory vector (maxcor=30), making the algorithm more inclined to utilize more distant gradient history information to achieve fine-grained local convergence. Both phases use the same L-BFGS-B algorithm core, differing only in the length of the memory vector.

[0190] After multiple iterations, the optimal solution vector for all MAP grid node parameters is output, which is the calibrated value (calibration value) of each MAP grid node. This optimization result can be directly written into the MAP table of the engine electronic control unit for actual vehicle control.

[0191] The engine MAP graph fitting optimization method provided in this application, in the current iteration step, determines gradient information based on gradient vector analysis, determines the search direction and step size based on boundary active set management and subspace minimization, updates the MAP grid node parameters according to the search direction and step size, and repeats the above steps until the convergence condition is met, thus obtaining the optimized MAP grid node parameters. The above method achieves stable and reliable solution to high-dimensional calibration problems under tight boundary constraints by using gradient vector analysis to drive iterative search and a projection gradient mechanism to ensure boundary feasibility, and obtaining the optimal MAP grid node parameters that satisfy the constraints under convergence conditions.

[0192] Figure 5 Example 4 of the engine MAP fitting optimization method provided for this solution, such as Figure 5 As shown, based on the foregoing embodiments, the method further includes:

[0193] S501: In the process of iterative optimization, the first round starts from the original MAP value, performs a complete iterative optimization, and records the obtained optimal solution as the safety net benchmark.

[0194] This step starts with the physically reasonable original MAP values ​​provided by the calibration engineer, performs the first round of complete iterative optimization, obtains an initial optimal solution, and uses it as a safety net benchmark to provide a reference anchor for comparison and backtracking in subsequent rounds, ensuring that the multi-starting point search process will not output a final solution that is worse than the result of the first round under any circumstances.

[0195] Specifically, the original MAP mesh node parameter values ​​provided by the calibration engineer are used as the initial starting point for the first round. These original MAP values ​​are derived from engineering experience or physical models and are inherently physically reasonable (i.e., the MAP surface is continuous, there is no jagged deformation, and the values ​​of each mesh point are within a reasonable physical range).

[0196] Starting with the original MAP value, execute a complete iterative optimization process in S203. This process includes: repeatedly calling S201-S202 to calculate the loss function value and analyze the gradient vector, iterating repeatedly through steps a to c of S203 (identifying the active set → determining the search direction in the free parameter subspace → determining the step size through line search → updating parameters and constraining boundary values ​​→ judging convergence) until the convergence condition is met, and obtaining the first round of optimal solution.

[0197] After the first round of iterative optimization, record the parameter vector corresponding to the optimal solution in the first round. and the loss function value corresponding to the solution This pair of data forms the safety net baseline, which is stored in memory for comparison and rollback in all subsequent rounds.

[0198] S502: In subsequent rounds, starting from the current optimal solution, the MAP grid is downsampled to a low resolution, random noise is applied to the low-resolution grid, and then the original resolution is restored by interpolation upsampling, generating multiple perturbation-based candidate starting points.

[0199] Starting from the current optimal solution, this step generates multiple spatially continuous perturbation candidate starting points through a process of downsampling, adding noise, and upsampling. This ensures that the perturbed MAP surface is naturally smooth and free from non-physical jagged deformation caused by independent noise, thereby exploring new parameter space regions under physically reasonable conditions and avoiding getting trapped in the current local optimum.

[0200] Specifically, the optimal solution obtained in the previous round (i.e., the current global optimal parameter vector) is used as the benchmark for the current round of perturbation, and multiple perturbation candidate starting points are generated based on this benchmark.

[0201] The current MAP grid node parameter matrix (original resolution, e.g., 20×20) is downsampled to obtain a low-resolution (e.g., 5×5) parameter matrix. Downsampling methods can employ average pooling or bilinear downsampling, merging parameter values ​​from adjacent regions in the high-resolution grid into a single value in the low-resolution grid.

[0202] Independent normal random noise is applied to each node of the low-resolution grid. The amplitude of this noise is controlled by the perturbation rate of the current round. In the early stages of optimization (such as the second round), the perturbation rate is set to a relatively high initial value (such as 20%), resulting in a larger perturbation amplitude, which is beneficial for exploration over a wider parameter space. As the optimization progresses and the error continues to improve, the perturbation rate is gradually reduced, resulting in a smaller perturbation amplitude, which is beneficial for fine-grained searching within the discovered promising regions.

[0203] The low-resolution parameter matrix after applying noise is restored to its original resolution (e.g., from 5×5 to 20×20) through interpolation upsampling. Bilinear interpolation can be used, ensuring that the perturbation of a noisy point in the low-resolution grid, after restoration to high resolution, smoothly and gradually affects its surrounding neighborhood. After upsampling, candidate starting points with perturbation are obtained.

[0204] Repeat the above steps multiple times, applying different random noise each time, to generate multiple perturbed candidate starting points (e.g., 3 candidate starting points) for subsequent screening.

[0205] S503: Select the one with the best error from multiple perturbation candidate starting points by performing forward propagation evaluation as the actual starting point, and re-execute a complete iterative optimization from the actual starting point.

[0206] This step involves rapidly evaluating multiple perturbation candidate starting points generated by S502 by performing only forward propagation (without performing complete iterative optimization and back propagation), selecting the one with the best error from the candidate starting points as the actual starting point for this round. This improves the quality of perturbation candidates with minimal computational cost. Then, complete iterative optimization is performed starting from this starting point, realizing a two-stage strategy of "rapid selection → fine optimization".

[0207] Specifically, the parameter vector of each candidate perturbation starting point is used as input, and forward propagation in S301 is performed once to calculate the loss function value corresponding to each candidate starting point (i.e., only the error between the predicted value and the measured value is calculated, and backpropagation is not performed to calculate the gradient). Since only forward propagation is performed, there is no need for backpropagation and iterative optimization. The computational cost of a single evaluation is much smaller than that of a complete iterative optimization (approximately milliseconds, which is negligible compared to the seconds to minutes of a complete optimization).

[0208] Compare the loss function values ​​of each candidate starting point, and select the candidate starting point with the smallest loss function value (i.e. the smallest fitting error) as the actual starting point for this round of optimization.

[0209] Using this actual starting point as the initial parameter, execute a complete iterative optimization process in S203 (i.e., call S201-S202 multiple times to calculate). The tuple is iterated through steps a-c until convergence to obtain the optimal solution for this round. After this round of optimization, the result is compared with the current global optimal solution (i.e., the safety net benchmark or the optimal solution in previous rounds). If the result of this round is better, the global optimal solution is updated.

[0210] S504: If the results obtained from subsequent iterations are worse than the safety net benchmark and the degree of deterioration exceeds the first preset threshold, then revert to the safety net benchmark and re-perturb.

[0211] When the optimization result of a certain round is worse than the safety net baseline and the degree of deterioration exceeds the first preset threshold, it indicates that the current perturbation direction may have caused the parameter space to enter a poor region, and continuing to search in this direction carries the risk of cascading failure. At this time, the safety net mechanism is triggered, and the perturbation is regenerated from the safety net baseline to prevent the multi-starting-point search strategy from outputting results worse than the first round due to continuous deterioration.

[0212] Specifically, after each round of optimization, the loss function value of the optimal solution obtained in this round is compared with the loss function value of the safety net benchmark:

[0213] Calculate the relative deterioration of the loss function value of the optimal solution in this round compared to the baseline loss function value of the safety net. Specifically, subtract the baseline loss function value of the safety net from the current round's loss function value, and then divide by the baseline loss function value of the safety net (or the absolute difference) to obtain the deterioration ratio.

[0214] If the degree of deterioration exceeds the first preset threshold, the result of this round is determined to be worse than the safety net benchmark and the deterioration is unacceptable, triggering the safety net rollback mechanism. If the degree of deterioration does not exceed the threshold (i.e., although the result of this round is slightly worse, it is still within the acceptable range), the safety net is not triggered, and the next round continues (perturbing from the optimal solution of this round).

[0215] When the safety net is triggered, all intermediate results from the current round and previous deterioration rounds are discarded, the parameters are reset to the parameter vector recorded by the safety net baseline, and a new perturbation candidate starting point is generated from this baseline (repeating S502), restarting the search process for subsequent rounds. This mechanism effectively prevents the cumulative performance degradation caused by continuously selecting inferior perturbation directions in multi-startpoint searches.

[0216] Optionally, during the iterative optimization process, if the error obtained from multiple consecutive optimization rounds does not decrease compared to the previous round, the perturbation rate is increased; if the error in the current round decreases by more than the second preset threshold compared to the previous round, the perturbation rate is decreased.

[0217] During the iterative optimization process, the perturbation rate is dynamically adjusted based on real-time feedback from the optimization status: when the error fails to decrease over multiple consecutive rounds, the perturbation rate is increased to escape the current local optimum; when the error continues to improve, the perturbation rate is decreased to refine the search within the currently promising region. This mechanism achieves automatic tuning of the search range throughout the entire process.

[0218] After each round of optimization, monitor the trend of the error in this round relative to the error in the previous round:

[0219] If the error obtained from multiple consecutive rounds of optimization (e.g., three consecutive rounds) does not decrease compared to the previous round, it indicates that the current search direction may have fallen into a local optimum and is difficult to improve further. In this case, increasing the perturbation rate (e.g., doubling it) will increase the noise amplitude in subsequent rounds, which will help to break out of the current local optimum and enter a wider parameter space for exploration.

[0220] If the reduction in error in this round compared to the previous round exceeds the second preset threshold (e.g., an improvement exceeding 5%), it indicates that the current search direction is correct and a promising descent path exists in the current basin. At this point, the perturbation rate is reduced (e.g., halved) to allow subsequent rounds to perform a more refined search within a smaller range, avoiding the destruction of the already found good parameter regions due to excessive perturbation.

[0221] Through the aforementioned adaptive mechanism, the disturbance rate is automatically adjusted throughout the process, with a large disturbance range in the exploration phase and a small disturbance range in the fine-tuning phase, thus achieving an automatic balance between "exploration" and "development".

[0222] Optionally, during the perturbation generation process of multi-starting point search, the perturbation intensity is differentiated according to the magnitude of the absolute gradient value of the parameters of each MAP grid node: parameters with large absolute gradient values ​​are subjected to smaller perturbation amplitudes, and parameters with small absolute gradient values ​​are subjected to larger perturbation amplitudes. This improves perturbation efficiency by focusing on exploring parameter regions that have not yet been fully optimized while maintaining the currently well-fitted region.

[0223] Specifically, during the optimization process, the absolute value of the gradient reflects the sensitivity of the loss function to changes in that parameter:

[0224] A large absolute gradient value means that even a small change in this parameter will cause a significant change in the loss function. In other words, this parameter is currently in a "sensitive" region, and its value plays a crucial role in reducing error. Applying a large perturbation to this parameter is likely to disrupt the already well-fitted state, causing a sharp rebound in error.

[0225] A small absolute gradient value means that a large change in this parameter has a limited impact on the loss function, indicating that the parameter is currently in an "insensitive" region and its value has not yet been fully optimized. Applying a larger perturbation at this point may allow us to explore a better value for the parameter without significantly damaging the already well-fitted region.

[0226] In practice, after each round of optimization, the analytical gradient vector g at the current optimal solution is obtained. This vector contains the gradient values ​​corresponding to all MAP grid node parameters. The absolute value of the gradient for each parameter is calculated, and parameters are sorted or proportionally allocated according to the magnitude of their absolute gradient values. For parameters with larger absolute gradient values, their perturbation amplitude coefficient is set to a smaller value (i.e., multiplied by a coefficient less than 1); for parameters with smaller absolute gradient values, their perturbation amplitude coefficient is set to a larger value (i.e., multiplied by a coefficient greater than 1). Finally, the perturbation amplitude of each parameter = base perturbation rate × perturbation amplitude coefficient corresponding to that parameter.

[0227] This mechanism complements the spatially correlated perturbation of downsampling-upsampling: downsampling-upsampling ensures the spatial continuity of the MAP surface after perturbation; gradient guidance ensures that the allocation of perturbation intensity is based on the current optimization state of the parameters, with high-sensitivity parameters protecting existing results and low-sensitivity parameters encouraging in-depth exploration. The synergistic effect of both mechanisms allows the multi-starting-point search strategy to further improve the quality and search efficiency of perturbation candidates while ensuring spatial smoothness.

[0228] The engine MAP fitting optimization method provided in this application, during the iterative optimization process, firstly starts from the original MAP value, performs a complete iterative optimization, and records the obtained optimal solution as a safety net benchmark. Subsequent rounds start from the current optimal solution, downsample the MAP grid to a low resolution, apply random noise to the low-resolution grid, and then upsample through interpolation to restore the original resolution, generating multiple perturbation candidate starting points. Through forward propagation evaluation, the candidate starting point with the best error is selected as the actual starting point. A complete iterative optimization is then performed again from the actual starting point. If the result of subsequent rounds of iterative optimization is worse than the safety net benchmark and the degree of deterioration exceeds a first preset threshold, the process reverts to the safety net benchmark and re-perturbs. This method, through safety net benchmark protection, spatially continuous perturbation generation, rapid forward propagation screening, and a safety net backoff mechanism, achieves physically reasonable multi-starting point search, improving the ability to escape local optima while avoiding cascading failures.

[0229] Figure 6 Example 5 of the engine MAP plot fitting optimization method provided for this solution, such as Figure 6 As shown, based on the foregoing embodiments, the method further includes:

[0230] S601: Receives user-inputted custom fitting expressions, parses the syntax structure of the expressions, extracts MAP lookup function calls, data column dependencies, and dependencies between Function nodes, and constructs a computational dependency topology in the form of a directed acyclic graph.

[0231] This step receives a custom fitting expression input by the user through a graphical interface. It extracts the composition structure and dependencies of the expression through lexical and syntactic parsing, and constructs a computational dependency topology in the form of a directed acyclic graph, providing a structured computational graph foundation for subsequent automatic differentiation and gradient calculation.

[0232] Specifically, users enter a custom fitting expression in the Function node editor of the graphical interface. This expression can include MAP lookup function calls (map_base(rpm,map_press)), data column references (such as data_column("torque_measured")), Function node references, and basic algebraic operations (addition, subtraction, multiplication, division, exponentiation, etc.). After receiving the expression string, it is preprocessed, including removing whitespace characters and identifying and marking various function calls, variable references, and operators.

[0233] Perform syntax parsing on the preprocessed expression string. The parsing process includes:

[0234] Lexical analysis: The expression string is broken down into a sequence of tokens, such as keywords (MAP function name, Function node name), identifiers (data column names, variable names), and operators (…). ), parentheses, and constant numbers, etc.

[0235] Syntax analysis: Based on the syntax rules of the expression, the token sequence is constructed into an abstract syntax tree (AST). Each tree node represents a unit of operation (such as a MAP lookup node, algebra operation node, function call node, data column reference node, etc.). The parent-child relationship between nodes reflects the hierarchical structure of the expression's operations.

[0236] Traverse the abstract syntax tree to extract three types of dependency relationships:

[0237] MAP lookup function calls: Identify all MAP lookup function calls in the expression, such as map_base, map_factor1, etc., extract the MAP name called and the corresponding input parameters (such as engine speed rpm, intake pressure map_press, intake temperature, ambient humidity, etc.), and record the position of each MAP lookup function node in the calculation graph and the input mapping relationship.

[0238] Data column dependencies: Identify references to measured data columns in the expression and establish the correspondence between the predicted value expression and the measured data;

[0239] Function node dependency relationships: Identify references to other Function nodes in expressions, extract the call hierarchy between Function nodes, determine which Function nodes are referenced by other nodes, which Function nodes exist independently, and whether there are circular dependencies, etc.

[0240] Based on the extracted dependencies, a computational dependency topology is constructed in the form of a directed acyclic graph (DAG):

[0241] Each node corresponds to a computational unit (such as MAP lookup function call node, algebraic operation node, Function node, data column reference node, loss function node, etc.).

[0242] Directed edges between nodes represent the direction of data dependency: the output of the predecessor node serves as the input of the successor node;

[0243] There are no directed cycles in the graph (all dependencies are unidirectional and there are no circular dependencies), ensuring that computation can proceed layer by layer from input to output in topological order;

[0244] The leaf nodes in the diagram are MAP grid node parameters (optimizable variables) and data column reference nodes (constant measured data), while the root node is the loss function node (final output).

[0245] This computational topology provides the structural foundation for subsequent dependency node identification and gradient propagation.

[0246] S602: Pre-compute non-dependent nodes in the computation-dependent topology, and retain expression strings for dependent nodes in the computation-dependent topology for subsequent dynamic evaluation.

[0247] This step is based on the computational dependency topology built by S601. Recursive dependency analysis is performed to identify Function nodes that the optimized target directly or indirectly depends on. Based on this, the nodes are divided into two categories: dependent nodes and non-dependent nodes. Different computational strategies are adopted for each category to balance computational efficiency and flexibility.

[0248] Specifically, recursive dependency analysis is performed on each Function node in the computational dependency topology: starting from the loss function node (root node), all its predecessor nodes are traversed in reverse along the edges of the computation graph; for each Function node encountered, it is checked whether the node appears on the reachable path starting from the loss function node; a breadth-first search (BFS) strategy is used to expand the search range layer by layer to identify all Function nodes that are directly or indirectly dependent on by the loss function node.

[0249] When a Function node is identified as a dependent node, all its input dependencies (i.e., the child Function nodes it references) are also marked as dependent nodes, ensuring the complete transitivity of dependencies.

[0250] Based on the results of recursive dependency analysis, the system categorizes nodes in the computational dependency topology into two classes:

[0251] Independent nodes: Function nodes that are not directly or indirectly dependent on by the loss function nodes. The computation results of these nodes do not participate in the calculation of the final loss function value, and therefore do not need to be recalculated in each iteration. A pre-computation is performed on these nodes during the preparation phase, and the computation results are cached as static tensors. In each subsequent forward propagation, this pre-computation result is directly read, without re-executing the computation logic.

[0252] Dependent nodes: Function nodes that are directly or indirectly dependent on by the loss function nodes. The computation results of these nodes directly affect the loss function value, so they need to participate in the computation and their computation path needs to be recorded in each forward propagation so that gradients can be traced during backpropagation. The expression strings of these nodes are preserved and dynamically evaluated in each forward propagation (i.e., the results are recalculated based on the current parameter values).

[0253] For Function nodes marked as dependent nodes, their original expression strings and associated parameter reference information are preserved. During each forward propagation, when the computation graph reaches this Function node, its expression string is dynamically evaluated based on the current input parameter values, and the result is output to its successor nodes. Since these nodes and their expressions fully participate in the forward propagation computation, their computation paths are recorded in the computation graph structure of the automatic differentiation framework, thus supporting gradient backpropagation during backpropagation.

[0254] For Function nodes marked as independent nodes, a pre-computation is performed during the preparation phase: using the data column reference values ​​required by the node expression as input (measured data constants), the output value of the Function node is calculated, and this output value is registered as a fixed constant (static tensor) in the computation graph. During the forward propagation of all subsequent iterations, this node no longer participates in computation, and its fixed output is directly used as the input of subsequent nodes.

[0255] S603: Replace the discontinuous branch selection function in the expression with a continuously differentiable transition function to obtain a computational graph of the traceable gradient propagation path.

[0256] This step replaces the discontinuous branch selection function (such as bswitch) used in the expression with a continuously differentiable sigmoid transition function (smooth_bswitch), ensuring that all nodes in the computation graph are differentiable, so that the automatic differentiation chain is not interrupted at the branch nodes, thereby ensuring that the analysis gradient vector can be calculated accurately.

[0257] Specifically, engine calibration expressions may contain conditional branching logic, such as composite formulas with saturation limiters (e.g., EGR_limit=bswitch(Y_1-30,Y_1,30), where bswitch determines the output Y_1 or 30 based on whether the intermediate variable Y_1 exceeds 30). In the abstract syntax tree of the expression, functions like bswitch are identified as discontinuous branch selection functions, whose output abruptly changes at the conditional threshold, the function value itself is discontinuous, and the derivative is either infinite or nonexistent.

[0258] Upon detecting a discontinuous branch selection function, it is automatically replaced with a continuously differentiable S-shaped transition function (smooth_bswitch). Specifically, the original bswitch(condition, true_value, false_value) outputs true_value when the condition is true and false_value when the condition is false, resulting in discontinuous jumps at the condition threshold.

[0259] The replaced smooth_bswitch uses an S-shaped transition curve (sigmoid function) to achieve a smooth transition between condition values, forming a continuous and differentiable transition band at the original discontinuous branch points.

[0260] The mathematical form of the S-shaped transition function guarantees that the function value is continuous on both sides of the conditional threshold and that there is a finite derivative at that point, so that the automatic differentiation can calculate the effective gradient value when passing through the node and continue to propagate the gradient back along the computation graph.

[0261] Through the above replacement, all nodes (including those that originally contained branching logic) become continuously differentiable functions. The entire computation graph, from input (MAP grid node parameters) to output (loss function value), forms a complete differentiable path. During backpropagation, the automatic differentiation framework can compute effective local gradients at each node, ensuring that the gradients propagate completely along the computation graph to all MAP grid node parameters without any interruptions.

[0262] After the dependency topology construction in S601, the node classification and branching strategy in S602, and the continuous replacement of the bswitch in S603, a complete computation graph with a traceable gradient propagation path is finally obtained. The core features of this computation graph include: recording the dependencies and topological order of all nodes; including references to pre-computed weight buffers for bilinear interpolation; all nodes are differentiable functions; gradients can propagate completely backward from the loss function node along the computation graph to all MAP grid node parameters; Function nodes marked as dependent nodes in S602 retain expression strings for dynamic evaluation, while nodes marked as non-dependent nodes store pre-computed constants. This computation graph will serve as the execution basis for forward propagation in S301 and backward propagation in S302, supporting gradient calculation throughout the optimization process. The pre-computed results of all non-dependent nodes and the expressions of dependent nodes in the computation graph are registered in the graph, ensuring direct access during each forward / backward propagation without repeated parsing or construction.

[0263] S604: For each MAP lookup function in the computation graph, pre-calculate the interpolation weights and indices for each data point, and register the interpolation weights and indices as fixed buffers for use in the lookup process during forward propagation.

[0264] This step pre-calculates the bilinear interpolation weights and one-dimensional flattened grid indexes for each data point in the computation graph based on its grid coordinate axes and the query coordinates of all data points. The weights and indexes are registered as fixed buffers, which degenerates the interpolation calculation in the forward propagation into pure table lookup multiplication and addition operations, significantly reducing the computational overhead of each iteration.

[0265] Specifically, for each MAP lookup function node in the computation graph, its grid coordinate axes (e.g., the X-axis is a list of engine speed values, and the Y-axis is a list of intake pressure values) are obtained from the definition information of that MAP. Then, the query coordinates of all data points are obtained from the dataset (i.e., the X-axis coordinate value and Y-axis coordinate value of each data point, such as the engine speed value and intake pressure value of each operating point).

[0266] For each data point's query coordinates (x, y), perform the following pre-calculation:

[0267] Locating the grid cell: Find the interval containing x in the X-axis coordinate list of the MAP. (Right now Find the interval containing y in the list of Y-axis coordinates. Determine the grid cell containing the data point, and the coordinates of the four corner points of that grid cell;

[0268] Calculate the interpolation weights: Calculate the interpolation coefficients w00, w10, w01, and w11 at the four corner points using the bilinear interpolation formula. Specifically, let... The weights of the four corner points are as follows:

[0269]

[0270]

[0271]

[0272]

[0273] Calculate the one-dimensional flattened grid index: convert the two-dimensional grid coordinates of the four corner points. Convert them to one-dimensional flat indexes respectively, the conversion method is as follows: ,get .

[0274] Allocate four weight buffers and four index buffers for each MAP lookup function node:

[0275] Weight buffer: Stores the interpolation weights w00, w10, w01, and w11 at the four corners of all data points. The size of each buffer is equal to the total number of data points, and the storage order is consistent with the order of the data points in the dataset.

[0276] Index Buffer: Stores the four corner points of all data points in a one-dimensional flat grid index idx00, idx10, idx01, and idx11. The size of each buffer is equal to the total number of data points, and the storage order is consistent with the order of data points in the dataset.

[0277] These buffers are created and filled all at once before optimization begins, and then registered to the corresponding MAP lookup function nodes in the computation graph, becoming fixed attributes of those nodes that remain unchanged throughout the optimization process without requiring recalculation. Simultaneously, the MAP grid node parameter values ​​are registered as optimizable parameters, allowing the automatic differentiation framework to track their gradients.

[0278] During the forward propagation of S301, when the calculation reaches the MAP lookup function node, coordinate positioning, interval search, or weight calculation are no longer performed. Instead, the weight and index values ​​corresponding to the current data point are directly read from the four weight buffers and four index buffers. The current parameter values ​​of the four corner points are extracted from the parameter vector of the current MAP grid node through four index lookups. The four corner point values ​​are multiplied by their corresponding weights and then summed to obtain the interpolation result for the data point. The entire interpolation process only involves four lookups, four multiplications, and three additions, without any search or dynamic calculation, significantly reducing the interpolation calculation time for a single forward propagation.

[0279] The engine MAP graph fitting optimization method provided in this application receives a user-input custom fitting expression, parses the expression's syntax structure, extracts MAP lookup function calls, data column dependencies, and dependencies between Function nodes, constructs a computational dependency topology in the form of a directed acyclic graph, pre-computes non-dependent nodes in the computational dependency topology, retains the expression strings for dependent nodes in the computational dependency topology for subsequent dynamic evaluation, replaces discontinuous branch selection functions in the expression with continuously differentiable transition functions, and obtains a computational graph with a traceable gradient propagation path. For each MAP lookup function in the computational graph, pre-computes the interpolation weights and indices for each data point, and registers the interpolation weights and indices as fixed buffers for lookup calls during forward propagation. This method, by parsing user-defined expressions to construct a differentiable computational graph, combined with node pre-computation branching and interpolation weight pre-computation buffering, achieves automatic differential calculation and efficient forward propagation of arbitrary fitting expressions, providing a traceable and complete path for accurate gradient calculation.

[0280] Optionally, at the start of each optimization run, by constructing a constraint configuration signature string, it is checked whether the optimal parameters saved in the previous run still meet the current boundary constraints. Only when both checks (constraint boundary out-of-bounds check + constraint signature consistency check) pass will the results of the previous run be inherited as the initial starting point; otherwise, it will restart from the original imported values. This maximizes the use of historical optimization results while ensuring the validity of parameters and prevents cross-run parameter contamination.

[0281] For example, a user obtains a better calibration result and saves the optimal parameters in an optimization run, but may subsequently modify the boundary constraints of the MAP (such as adjusting the upper and lower limits of a grid point, changing the configuration of a function node, etc.). If the previously saved optimal parameters are directly inherited as the starting point for the next optimization without any checks, these parameters may no longer meet the new constraints, causing the optimization process to start from an "illegal" starting point, leading to optimization instability or unusable results.

[0282] Specifically, at the start of the optimization run, all constraint configuration information for the current optimization task is collected, including but not limited to: the upper and lower bounds of all MAP grid node parameters (lo and hi arrays); the upper and lower bounds of the user-defined gradient constraints (maximum / minimum gradient values ​​in the X and Y directions); the configuration parameters of the Function node; the lock points and their lock values, etc.

[0283] These constraint configuration information are concatenated into a string according to a fixed field order and format. Then, a hash calculation (such as MD5 or SHA-1) is performed on this string to generate a fixed-length constraint configuration signature (hereinafter referred to as "constraint signature"). The core characteristic of this signature is that the same constraint configuration will inevitably produce the same signature, and any modification to the constraint configuration (regardless of the magnitude) will cause a significant change in the signature.

[0284] At the end of each optimization run, the optimal parameter vector and corresponding constraint signatures of this run are saved to persistent storage (such as a local file or database). The saved optimal parameters are the historical best-seen parameters tracked internally by the optimizer, rather than potentially worse final parameters returned by the optimizer due to a last-step bounce or numerical issues, ensuring that the saved parameters are always the optimal solutions obtained in this run.

[0285] At the start of each optimization run, perform the following two checks:

[0286] The first level of detection—boundary violation check: Read the optimal parameter vector saved from the last run and check each parameter to see if it is still within the current boundary range. If any parameter exceeds the current boundary, the first detection fails, and the system determines that the previous result cannot be inherited.

[0287] The second layer of checks—constraint signature consistency checks—constructs the constraint configuration signature for the current run and compares it with the constraint signature saved from the previous run. If the two signatures are completely identical, the second layer of checks passes; if the signatures are inconsistent (indicating that the constraint configuration has been modified), the second layer of checks fails, and the previous result is determined not to be inherited.

[0288] Only when both checks pass are the optimal parameters from the previous run inherited as the initial starting point for this optimization, and the optimizer continues to iterate from that starting point. If either check fails, the previously saved optimal parameters are discarded, and the optimization is restarted from the original imported values ​​provided by the calibration engineer, ensuring that each optimization starts from a physically reasonable starting point that satisfies all current constraints.

[0289] This mechanism effectively prevents the following scenarios of pollution across operating parameters:

[0290] For example, if a user modifies the lower bound of a MAP grid point (e.g., from 0 to 5), but the previously saved parameter for that point was 3, directly inheriting it will result in an invalid parameter.

[0291] The user modified the gradient constraints (e.g., changed the maximum rate of change from 10 to 5), but there are adjacent points with a difference value of 8 in the previously saved parameters. Directly inheriting these parameters will cause subsequent optimizations to produce intermediate results that violate the new constraints.

[0292] Users run optimizations under different constraint configurations, and the saved optimal parameters are mixed and used interchangeably.

[0293] After optimization, multi-panel, multi-dimensional analysis and evaluation tools are provided to comprehensively diagnose the fit quality and MAP characteristics, helping calibration engineers to examine the optimization results layer by layer from macro statistics to micro residuals, and complete the fit quality assessment in one stop without writing additional analysis code.

[0294] Specifically, based on a unified vectorized prediction calculation engine, a complete prediction calculation is performed on all data points after optimization (using the optimized MAP grid node parameters, bilinear interpolation is performed on the query coordinates of each data point to obtain the prediction value). This prediction calculation result serves as the basic data source shared by all analysis panels, avoiding independent and repetitive calculations by each panel and improving analysis response speed.

[0295] Design five analysis panels:

[0296] (1) Statistical summary panel

[0297] This panel displays key statistical metrics for fit quality, including: Root Mean Square Error (RMSE) and Mean Absolute Error (MAE), used to assess overall fit accuracy; Maximum Absolute Error and error percentile distributions (e.g., P50, P90, P95, P99), used to assess the central tendency and long tail of the error distribution; and the coefficient of determination (…). The MAP model is used to evaluate how well it interprets the variance of the measured data; it provides runtime information such as the number of valid data points, the number of MAPs involved in the fitting, the total number of iterations, and the total computation time.

[0298] (2) Data comparison panel

[0299] The panel lists the measured value, predicted value, and residual (measured value - predicted value) for each data point in a table format, row by row. It supports sorting by the absolute value of the residual (from largest to smallest or smallest to largest), which helps engineers quickly locate the data point with the largest deviation and check whether there is a measurement anomaly or whether the data point is located in the MAP extrapolation area.

[0300] (3) MAP range view panel (including heat map)

[0301] This panel visualizes three sets of MAP surfaces in the form of heatmaps: Original surface: the initial MAP surface provided by the calibration engineer (parameter value distribution before optimization); Fitted surface: the optimal MAP surface obtained after optimization (parameter value distribution after optimization), intuitively showing the effect of optimization on the MAP surface; Residual surface: a heatmap of the difference between the fitted surface and the original surface (fitted value - original value), with color intensity indicating the adjustment magnitude and direction of parameter values ​​at each grid point, helping engineers identify which areas were adjusted the most by the optimizer and determine whether the adjustment meets physical expectations.

[0302] Each view provides interactive features (such as zooming, panning, and adjusting color levels).

[0303] (4) Scatter verification plot panel

[0304] This panel displays the correspondence between predicted values ​​(Y-axis) and measured values ​​(X-axis) in scatter plot format, with each data point corresponding to a scatter plot marker. Ideally, all points should fall near the 45-degree diagonal (y=x). This panel also provides:

[0305] The fitted regression line (linear regression) is used to assess the systematic bias between the prediction and the actual result (whether the slope of the regression line is close to 1 and whether the intercept is close to 0).

[0306] The R² and RMSE values ​​are marked in a prominent position on the graph.

[0307] (5) Residual distribution analysis panel

[0308] This panel provides statistical and visual analysis of the residual data, including:

[0309] The residual histogram and the superimposed normal distribution fitting curve are used to check whether the residuals approximately follow a normal distribution (if the residual distribution deviates significantly from the normal distribution, it may indicate that there is a systematic bias in the fitting model).

[0310] A scatter plot of residuals versus measured values ​​is used to check whether the residuals are evenly distributed on both sides of the zero line and whether they show a trend change with the size of the measured values ​​(if the residuals increase with the increase of the measured values, it suggests that there may be heteroscedasticity).

[0311] Scatter plots of residuals versus input parameters are used to check whether the residuals exhibit clustered biases in specific regions of the input parameter space, helping engineers identify problems with poor MAP fitting in certain operating conditions.

[0312] Using the five analysis panels mentioned above, engineers can quickly and comprehensively evaluate the quality and reliability of the optimization results and determine whether it is necessary to adjust the objective function weights, constraint boundaries, or expression structure before re-optimizing.

[0313] For example, the method provided in this application is illustrated using the MAP map fitting optimization of an 825-parameter engine as an example. This embodiment involves the joint optimization of three MAP parameters: map_base (20 rows × 20 columns = 400 parameters), map_factor1 (20 rows × 20 columns = 400 parameters), and map_factor2 (5 rows × 5 columns = 25 parameters), totaling 825 optimizable parameters. The fitting formula is a typical cascaded correction structure:

[0314]

[0315] Step 1: Importing and Standardizing Multi-Source Data

[0316] The data import dialog box on the data import page allows users to select data files in various formats, including standard description file format for engine electronic control units, tabular data format (CSV), spreadsheet format (Excel), and plain text format (TXT).

[0317] For the standard description file format of the engine electronic control unit, its structured content is parsed: reading the data block describing the two-dimensional grid size, reading the axis vector data block, and reading the grid point value matrix block. Data encoding is automatically identified during the parsing process to ensure correct display of language tags.

[0318] For tables / spreadsheets / plain text formats, a configuration dialog box allows users to specify parameters such as delimiters, header rows, and encoding, loading table data into a unified data structure.

[0319] Each imported MAP data is created as a MAP parameter object, whose attributes include: row axis vectors (e.g., 20 breakpoint values ​​for engine speed RPM), column axis vectors (e.g., 20 breakpoint values ​​for load), grid point value matrix (a 20×20 array), a deep copy of the original imported values ​​(as a reference value for cross-operation protection), a set of lockpoint coordinates, and a mapping dictionary from lockpoint coordinates to lockpoint values.

[0320] For the imported training data (measured data points), it is loaded as a dataset object. Its data structure contains multiple columns of measurement data, with typical columns including: rotational speed, load, and measured output value.

[0321] Step 2: Editing the user-defined fitting formula

[0322] Users create FunctionNode objects in the Function node editing page. Each FunctionNode contains an expression string that supports the following mathematical operators and functions: arithmetic operations, exponentiation, parenthesis precedence control, square root, absolute value, exponentiation, logarithm, sine, cosine, tangent, power function, mean, summation, and smooth differentiable binary switching functions.

[0323] Users use a specific syntax when referencing MAP parameters in expressions. During expression parsing, these MAP references are automatically identified and replaced with bilinear interpolation calculation code.

[0324] Users configure optimization parameters for each MAP parameter on the parameter configuration page: optimization switch, lower / upper bound, smoothness weight S, and gradient constraints (such as upper and lower bounds of gradient in the x-direction and gradient in the y-direction).

[0325] Users can set optimization hyperparameters on the optimization page: number of iterations, number of rounds with multiple starting points, and tolerance.

[0326] Step 3: Pre-calculation of Bilinear Interpolation Weights and Indices

[0327] Extract the axis coordinates of each data point from the training data.

[0328] For each MAP parameter, perform the following one-time pre-calculation for each data point:

[0329] Locate the interval of the data point in the x-axis vector based on its x-coordinate, and record the left index and the normalized position within the interval; locate the interval of the data point in the y-axis vector based on its y-coordinate, and record the index and the normalized position within the interval; calculate the bilinear interpolation weights of the four corner points; calculate the four one-dimensional flattened grid indices.

[0330] Weights and indices are stored in a fixed buffer and are automatically managed by the model; they do not participate in gradient calculation. This pre-computation is performed only once.

[0331] Step 4: Forward Propagation Calculation

[0332] Construct a differentiable objective function encapsulation model. Its core data members include: a MAP parameter dictionary, a pre-computed weight buffer, and a pre-computed index buffer.

[0333] For each data point, bilinear interpolation is performed on the grid point values ​​of each MAP parameter. The parameter grid value matrix is ​​flattened into a one-dimensional vector. The four grid cell values ​​are retrieved from the one-dimensional vector using a pre-computed flattening index, multiplied by four pre-computed weight coefficients, and then summed. A single interpolation calculation involves four table lookups, four multiplications, and three additions, with zero search and zero conditional branches.

[0334] The bilinear interpolation results of all MAP parameters are substituted as local variables into user-defined expressions for computation. The expressions are executed in a protected computational environment, with built-in function mappings to automatic differentiating operators, ensuring all computations are performed on a differentiable computational graph. Specifically, discontinuous branch selection functions are mapped to sigmoid smooth transition functions with minimal forward error and non-zero gradients maintained in the transition region.

[0335] The expression yields the predicted value vector, which is then compared with the measured value vector to calculate the sum of squared errors and the root mean square error.

[0336] For each MAP parameter involved in the optimization, a smoothing penalty term is calculated. For a shape of... The MAP parameter mesh is used to calculate the second-order central difference and its cross terms in the x and y directions, making the penalty stronger when the MAP surface is bowl-shaped or dome-shaped, effectively suppressing local high-frequency problems caused by bending in the same direction.

[0337] For each MAP parameter containing a lockpoint, calculate the concentric circle transition penalty for the lockpoint. For each lockpoint, define three concentric rings using Manhattan distance as the metric: the weights of the four direct neighbors in the first ring (Manhattan distance = 1) are 1.0, the weights of the eight neighbors in the second ring (Manhattan distance = 2) are 0.5, and the weights of the twelve neighbors in the third ring (Manhattan distance = 3) are 0.25. For each neighbor on a ring, calculate the penalty term: total lock penalty = smoothness weight × fixed coefficient × sum of penalties from all rings.

[0338] For each MAP parameter configured with gradient constraints, calculate the gradient constraint penalty. For the first-order difference in the x and y directions, check if it exceeds the user-defined upper and lower bounds of the gradient. Apply a secondary penalty to the portion exceeding the bounds, and apply no penalty to the portion within the bounds.

[0339] The total loss function is the sum of the terms: sum of squared errors + smoothing penalty + lockpoint concentric circle transition penalty + gradient constraint penalty.

[0340] Step 5: Accurately analyze gradient calculation

[0341] After forward propagation is complete, the model's computation graph contains the complete computational chain from 825 parameters to the total loss value.

[0342] The backpropagation interface of the automatic differentiation framework is called, and the computation graph is traversed through backpropagation. The chain rule is used to accurately calculate the partial derivative of each parameter with respect to the total loss, resulting in an 825-dimensional analytical gradient vector.

[0343] Return a tuple consisting of the gradient vector and the loss function value.

[0344] Step Six: Projection Gradient Boundary Constraint Optimization

[0345] Define a module that encapsulates the objective function. This module executes the following process each time it is called by the optimizer:

[0346] (1) Convert the input NumPy parameter vector into a tensor format of the automatic differentiation framework;

[0347] (2) Execute safety clamping to limit the parameters within the preset boundary range;

[0348] (3) Call the model parameter loading interface to load the parameters into the model (and trigger the lock point forced recovery: traverse all locked mesh cells and force their values ​​to be restored to the corresponding locked values).

[0349] (4) Use forward propagation to obtain the loss function value and root mean square error;

[0350] (5) Use backpropagation to obtain the analysis gradient vector;

[0351] (6) Internal tracking of the best root mean square error in history and the corresponding parameter vector;

[0352] (7) Returns the loss function value and gradient vector tuple.

[0353] Call the interface of the quasi-Newton optimizer with boundary constraints:

[0354] Objective function: Return tuple;

[0355] Initial parameter vector: 825 dimensions;

[0356] Optimization method: Quasi-Newton method with boundary constraints;

[0357] Gradient mode: Objective function returns tuple;

[0358] Boundary list: the upper and lower bounds of each element;

[0359] Options include: L-BFGS-B memory vector logarithm, function evaluation budget, number of iterations, function value convergence tolerance, and gradient convergence tolerance.

[0360] The optimizer's internal execution process:

[0361] (1) Using a two-loop recursive algorithm to approximate the inverse of the Hessian matrix, only the nearest m pairs of vectors need to be stored, resulting in a space complexity of O(825×m) instead of ;

[0362] (2) Automatically check boundary conditions in each line search and perform projection operation on out-of-bounds trial points (the boundary constraint mechanism built into the optimizer).

[0363] (3) Terminate when the function evaluation budget is exhausted or the convergence tolerance condition is met.

[0364] Post-processing is performed after optimization:

[0365] (1) Securely clamp the parameters returned by the optimizer;

[0366] (2) Calculate the final root mean square error;

[0367] (3) Compare the root mean square error of the historical best with the root mean square error of the optimizer at the last iteration point. If the historical best is better, then the historical best parameters will be used to cover the final result. This optimal parameter protection mechanism effectively defends against error bounce caused by the failure of line search or the step size being too small in the last step of the quasi-Newton method.

[0368] Anomaly recovery: When the optimizer throws an exception due to numerical issues, the exception is captured through the exception capture mechanism, and the system automatically rolls back to the historical best parameters, ensuring that the optimization process does not lose all the progress already made due to a single crash.

[0369] Step 7: Multi-starting point spatial correlation perturbation strategy

[0370] (Round 1): The first round starts from the original parameter values ​​imported by the user and performs a complete quasi-Newton method optimization with boundary constraints. After optimization, the optimal solution vector of the first round is saved as the safety net baseline, and the global optimal solution is initialized.

[0371] Round 2: Initialize the perturbation rate to 20%;

[0372] Round 3: The perturbation rate is initialized to 12%;

[0373] Third round and above: Dynamically adjusted based on historical improvement;

[0374] If the root mean square error does not improve after several consecutive rounds: double the perturbation rate and increase the search range;

[0375] If the improvement exceeds the preset threshold: the disturbance rate is halved, and the search range is narrowed to finely optimize the current basin.

[0376] Gradient-guided perturbation intensity allocation, applicable to rounds three and beyond:

[0377] (1) Perform an additional forward and backward propagation at the global optimal solution to obtain an 825-dimensional gradient vector;

[0378] (2) Calculate the absolute value of the gradient;

[0379] (3) Take a high quantile value of the absolute value of the gradient as a reference scale;

[0380] (4) Calculate the disturbance intensity scaling factor for each parameter;

[0381] (5) Meaning: Parameters with large absolute gradient values ​​(the loss is highly sensitive to these parameters) receive smaller perturbation amplitudes, avoiding the destruction of existing good fits; parameters with small absolute gradient values ​​receive larger perturbation amplitudes, encouraging exploration.

[0382] Spatial-related disturbance generation:

[0383] (1) For shapes that are Use the MAP parameters to create the downsampling grid size;

[0384] (2) Generate independent normally distributed random noise on the downsampling grid;

[0385] (3) Upsample the noise to full resolution using interpolation;

[0386] (4) Full-resolution noise multiplied by noise standard deviation (determined by perturbation rate and parameter boundary range) and gradient-guided scaling factor;

[0387] (5) Generate perturbation candidate parameters and clamp them to the boundary range.

[0388] The key advantage of spatially correlated perturbation is that the downsampling-upsampling process ensures that the perturbed MAP surface maintains spatial continuity and is naturally smooth, avoiding non-physical jagged deformation that may occur when noise is added independently grid by grid on the full-resolution grid.

[0389] Multiple candidate quick filtering:

[0390] (1) For each perturbation candidate, a forward propagation-only evaluation is performed, calculating only the root mean square error and not the gradient. Each evaluation takes a very short time.

[0391] (2) Generate multiple perturbation candidates in each round and evaluate their root mean square error values ​​respectively;

[0392] (3) Set the root mean square error threshold, that is, the error of the perturbation candidate must be lower than a certain multiple of the current best error and not exceed a certain absolute upper limit;

[0393] (4) If the error of the best candidate among multiple candidates is lower than the threshold, then select the candidate as the starting point for the next round of optimization;

[0394] (5) If all candidates exceed the threshold, the perturbation rate will be automatically halved and candidates will be regenerated, with a maximum of several retries. Even if the threshold is still exceeded after the number of retries is exhausted, the relatively optimal candidate will be selected to avoid the predicament of having no candidates available.

[0395] Safety net mechanism:

[0396] When multiple consecutive rounds of optimization results are deemed invalid (i.e., the root mean square error after optimization is worse than the global optimum before perturbation, and the degree of deterioration exceeds a preset threshold), a safety net is triggered: the current global optimum is reset to the result of the first round, and the next round's starting point is generated by re-perturbing from the first round's result. This mechanism prevents the "cascading bad luck" that may result from multiple rounds of perturbation, i.e., the situation where continuous perturbation from a suboptimal solution leads to getting further and further away from the true optimum.

[0397] Step 8: In-depth fine-tuning

[0398] After all multi-starting point rounds have been executed, select the parameter vector corresponding to the global optimal solution.

[0399] The optimization of the global optimal solution is performed in a deep fine-tuning stage, with different parameter settings:

[0400] Larger L-BFGS-B memory vectors: better approximate the Hessian inverse matrix structure of the 825-dimensional parameter space;

[0401] More function evaluation budget: allows for more complete convergence;

[0402] More stringent convergence tolerance.

[0403] The purpose of deep fine-tuning is to achieve higher accuracy in local optima within the identified high-quality basins. The effectiveness of the two-stage architecture lies in the fact that multiple starting points are responsible for coarse-grained global search, while deep fine-tuning is responsible for achieving high-accuracy local convergence within the optimal basins using a larger evaluation budget and stricter tolerance.

[0404] Record the initial root mean square error before deep fine-tuning. If the root mean square error after deep fine-tuning deteriorates by more than a preset proportion compared to the initial value, revert to the parameter value before deep fine-tuning.

[0405] Step Nine: Post-processing and Diagnostic Analysis

[0406] The final optimal parameter vector is clamped element by element to the boundary range, eliminating minor out-of-bounds errors that may occur during the iteration process due to numerical errors.

[0407] Iterate through all lock points of all MAP parameters and force their values ​​to be restored to the original lock values, correcting even the smallest deviations.

[0408] Perform a complete forward propagation using the final parameters and calculate the root mean square error (RMSE) as the core indicator of the fitting quality.

[0409] Set the smoothness weight to 0 and recalculate the root mean square error. If the zero-smoothness root mean square error is almost the same as the normal root mean square error, it indicates that the smoothing penalty has not had a substantial impact on the optimization, suggesting that users can try reducing the smoothness weight value to obtain a lower root mean square error.

[0410] The built-in multi-panel statistical diagnostic module automatically generates the following analysis results:

[0411] Panel 1 - Statistical Summary: Displays statistical indicators of measured and predicted data side-by-side: minimum, maximum, mean, standard deviation, range, and number of data points.

[0412] Panel 2 - Data Table: A complete list of data points, including the measured value, predicted value, residual, and percentage of relative error for each data point.

[0413] Panel 3 - Range View: After the user selects the MAP parameter to view, three two-dimensional heatmaps are displayed side by side: the original imported MAP value, the optimized MAP value, and the difference between the optimized and original values. The heatmaps use a uniform color mapping scale.

[0414] Panel 4 - Scatter Plot: A scatter plot of measured and predicted values, with the horizontal axis representing measured values ​​and the vertical axis representing predicted values. An ideal reference line and a linear regression line are overlaid to provide numerical annotations for the regression equation and coefficient of determination.

[0415] Panel 5 - Residual Analysis: Includes three subplots: residual distribution histogram and normal distribution fitting curve, residual and predicted value scatter plot (heteroscedasticity detection), and residual and data sequence scatter plot (trend detection).

[0416] All panels share a unified vectorized prediction engine, which calculates the predicted values ​​of all data points in batches at once based on pre-calculated bilinear interpolation weights and recursive node expressions. Multiple panels can reuse the same prediction calculation results to avoid redundant calculations.

[0417] Figure 7 A schematic diagram of the engine MAP fitting and optimization device provided for this solution is shown below. Figure 7 As shown, the engine MAP fitting and optimization device 70 specifically includes:

[0418] Function construction module 71 is used to construct an objective function, which includes an error term and at least one penalty term for constraining the characteristics of the MAP parameters;

[0419] The computation module 72 is used to perform one forward propagation and one backward propagation on the objective function based on a pre-built computation graph through an automatic differentiation framework, to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters;

[0420] The iterative optimization module 73 is used to iteratively optimize the MAP grid node parameters based on the loss function value and the analytical gradient vector using a boundary-constrained optimization algorithm based on the analytical gradient. In each iteration, the parameter components that touch the boundary are constrained to the boundary value by the boundary active set management, and the parameter update is performed in the free parameter subspace to obtain the optimized MAP grid node parameters.

[0421] Furthermore, the boundary-constrained optimization algorithm based on analytical gradients is a boundary-constrained quasi-Newton optimization algorithm; performing parameter updates within the free parameter subspace is equivalent to performing quasi-Newton optimization within the free parameter subspace.

[0422] Furthermore, the calculation module 72 is specifically used for:

[0423] Using the MAP grid node parameters as input, forward propagation is performed along the computation graph to obtain the loss function value;

[0424] Perform a backpropagation along the computation graph, and use the chain rule to propagate the gradient of the loss function value back layer by layer to obtain the analysis gradient vector.

[0425] Furthermore, the penalty items include a concentric circle transition penalty item for the lock point;

[0426] The concentric circle transition penalty term for the lock point divides the neighborhood into multiple rings based on the Manhattan distance from each lock point. Neighborhood points in different rings are assigned different attenuation weights, and the attenuation weights decrease from the inside out in each ring. This results in a greater penalty for neighboring points that are closer to the lock point. This is used to apply smooth transition constraints to the free parameters around the user-selected lock grid point.

[0427] Furthermore, the engine MAP fitting and optimization device 70 also includes:

[0428] The hard constraint module 74 is used to apply an immovable hard constraint to the locked grid points by setting the upper and lower bounds of the parameter elements corresponding to the locked grid points to the tolerance range of the locked values ​​when constructing the parameter boundaries of the objective function, in response to the locked grid points and locked values ​​selected by the user on the MAP diagram.

[0429] Furthermore, the engine MAP fitting and optimization device 70 also includes:

[0430] The forced recovery module 75 is used to forcibly restore the parameter values ​​corresponding to the user-selected locked grid points to the locked values ​​after each parameter update and when the final result is output during the iterative optimization of MAP grid node parameters.

[0431] Furthermore, the penalty terms include smoothing penalty terms;

[0432] The smoothing penalty term consists of the weighted sum of squared second-order differences in the X direction, the weighted sum of squared second-order differences in the Y direction, and cross terms, and is used to constrain the smoothness of the MAP surface.

[0433] Furthermore, the penalty term includes a gradient constraint penalty term;

[0434] The gradient constraint penalty term applies a unidirectional penalty to the first-order difference that exceeds the preset upper and lower limits of the gradient, which is used to constrain the rate of change of MAP in the coordinate axis direction.

[0435] Furthermore, the iterative optimization module 73 specifically includes:

[0436] Step a: In the current iteration step, determine the gradient direction based on the analyzed gradient vector;

[0437] Step b: Based on the gradient direction, in each iteration, the parameter components that touch the boundary are constrained to the boundary values ​​by the boundary active set management, and quasi-Newton optimization is performed in the free parameter subspace to determine the search direction and step size.

[0438] Step c: Update the MAP grid node parameters according to the search direction and step size;

[0439] Repeat steps a-c above until the convergence condition is met to obtain the optimized MAP grid node parameters.

[0440] Furthermore, the engine MAP fitting and optimization device 70 also includes:

[0441] Multi-starting point search module 76, used for:

[0442] During the iterative optimization process, the first round starts from the original MAP value, performs a complete iterative optimization, and records the obtained optimal solution as the safety net benchmark.

[0443] In subsequent rounds, starting from the current optimal solution, the MAP grid is downsampled to a low resolution, random noise is applied to the low-resolution grid, and then the original resolution is restored through interpolation upsampling, generating multiple perturbation-based candidate starting points;

[0444] By performing forward propagation evaluation, the error-optimized one from multiple perturbation candidate starting points is selected as the actual starting point, and a complete iterative optimization is performed again from the actual starting point;

[0445] If the results obtained from subsequent iterations are worse than the safety net benchmark and the degree of deterioration exceeds the first preset threshold, then revert to the safety net benchmark and re-perturb.

[0446] Furthermore, the multi-starting-point search module 76 is also used for:

[0447] If the error obtained from multiple consecutive optimization rounds does not decrease compared to the previous round, then increase the perturbation rate;

[0448] If the error in this round decreases by more than the second preset threshold compared to the previous round, then the perturbation rate is reduced.

[0449] Furthermore, the engine MAP fitting and optimization device 70 also includes:

[0450] Graph building module 77, used for:

[0451] It receives user-inputted custom fitting expressions, parses the syntax structure of the expressions, extracts MAP lookup function calls, data column dependencies, and dependencies between Function nodes, and constructs a computational dependency topology in the form of a directed acyclic graph.

[0452] Pre-compute non-dependent nodes in the computation-dependent topology, and retain expression strings for dependent nodes in the computation-dependent topology for subsequent dynamic evaluation;

[0453] By replacing the discontinuous branch selection function in the expression with a continuously differentiable transition function, a computational graph of the traceable gradient propagation path is obtained.

[0454] Furthermore, the engine MAP fitting and optimization device 70 also includes:

[0455] Interpolation module 78 is used to calculate each MAP lookup function in the graph, pre-calculate the interpolation weights and indices for each data point, and register the interpolation weights and indices as fixed buffers for lookup calls during forward propagation.

[0456] Figure 8 The schematic diagram of the electronic device structure provided for this solution is as follows: Figure 8 As shown, the electronic device 80 includes a memory 81 and a processor 82;

[0457] Memory 81 stores instructions executed by the computer;

[0458] The processor 82 executes computer execution instructions stored in the memory 81 to implement the methods provided in the above embodiments.

[0459] Optionally, the electronic device 80 also includes a communication component 83. The processor 82 and the communication component 83 are connected via a bus 84.

[0460] In a specific implementation, at least one processor 82 executes computer execution instructions stored in memory 81, causing at least one processor 82 to perform the above-described method.

[0461] The specific implementation process of processor 82 can be found in the above method embodiments, and its implementation principle and technical effect are similar. It will not be repeated here.

[0462] In the above embodiments, it should be understood that the processor can be a Central Processing Unit (CPU), or other general-purpose processors, digital signal processors (DSPs), application-specific integrated circuits (ASICs), etc. The general-purpose processor can be a microprocessor or any conventional processor. The steps of the method disclosed in this invention can be directly implemented by a hardware processor, or implemented by a combination of hardware and software modules within the processor.

[0463] The memory may include random access memory (RAM) and may also include non-volatile memory (NVM), such as at least one disk storage device.

[0464] The bus can be an Industry Standard Architecture (ISA) bus, a Peripheral Component Interconnect (PCI) bus, or an Extended Industry Standard Architecture (EISA) bus, etc. Buses can be categorized as address buses, data buses, control buses, etc. For ease of illustration, the buses shown in the accompanying drawings are not limited to a single bus or a single type of bus.

[0465] This application also provides a computer program product, including a computer program that, when executed by a processor, implements the above-described method.

[0466] This application also provides a computer-readable storage medium storing computer-executable instructions, which, when executed by a processor, implement the above-described method.

[0467] The aforementioned readable storage medium can be implemented by any type of volatile or non-volatile storage device or a combination thereof, such as static random access memory (SRAM), electrically erasable programmable read-only memory (EEPROM), erasable programmable read-only memory (EPROM), programmable read-only memory (PROM), read-only memory (ROM), magnetic storage, flash memory, magnetic disk, or optical disk. The readable storage medium can be any available medium accessible to a general-purpose or special-purpose computer.

[0468] An exemplary readable storage medium is coupled to a processor, enabling the processor to read information from and write information to the readable storage medium. Of course, the readable storage medium can also be a component of the processor. The processor and the readable storage medium can reside in an Application Specific Integrated Circuit (ASIC). Alternatively, the processor and the readable storage medium can exist as discrete components in the device.

[0469] The division of units is merely a logical functional division; in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the coupling or direct coupling or communication connection shown or discussed may be indirect coupling or communication connection through some interfaces, devices, or units, and may be electrical, mechanical, or other forms.

[0470] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0471] In addition, the functional units in the various embodiments of the present invention can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0472] If a function is implemented as a software functional unit and sold or used as an independent product, it can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of this invention, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods of the various embodiments of this invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0473] Those skilled in the art will understand that all or part of the steps of the above-described method embodiments can be implemented by hardware related to program instructions. The aforementioned program can be stored in a computer-readable storage medium. When executed, the program performs the steps of the above-described method embodiments; and the aforementioned storage medium includes various media capable of storing program code, such as ROM, RAM, magnetic disks, or optical disks.

[0474] The above embodiments are merely preferred embodiments provided to fully illustrate the present invention, and the scope of protection of the present invention is not limited thereto. Equivalent substitutions or modifications made by those skilled in the art based on the present invention are all within the scope of protection of the present invention.

Claims

1. A method for fitting and optimizing engine MAP diagrams, characterized in that, include: Construct an objective function, which includes an error term and at least one penalty term for constraining the characteristics of the MAP parameters; Based on a pre-built computational graph, the objective function is subjected to one forward propagation and one backward propagation through an automatic differentiation framework to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters. Based on the loss function value and the analysis gradient vector, the MAP grid node parameters are iteratively optimized using an analysis gradient-based bounded constraint optimization algorithm. In each iteration, the parameter components that touch the boundary are constrained to the boundary value by the boundary active set management, and the parameter update is performed in the free parameter subspace to obtain the optimized MAP grid node parameters.

2. The method according to claim 1, characterized in that, The boundary-constrained optimization algorithm based on analytical gradients is a boundary-constrained quasi-Newton optimization algorithm; the parameter update in the free parameter subspace is equivalent to performing quasi-Newton optimization in the free parameter subspace.

3. The method according to claim 1, characterized in that, Based on the pre-constructed computational graph, the objective function is subjected to one forward propagation and one backward propagation using an automatic differentiation framework to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters, including: Using the MAP grid node parameters as input, forward propagation is performed along the computation graph to obtain the loss function value; Perform a backpropagation along the computation graph, and use the chain rule to propagate the gradient of the loss function value back layer by layer to obtain the analysis gradient vector.

4. The method according to claim 1, characterized in that, The penalty items include a concentric circle transition penalty item at the lock point; The concentric circle transition penalty term for the locking point divides the neighborhood into multiple rings based on the Manhattan distance from each locking point. Neighborhood points in different rings are assigned different attenuation weights, and the attenuation weights decrease from the inside out in each ring, so that the closer the neighboring point is to the locking point, the greater the penalty it receives. This is used to apply smooth transition constraints to the free parameters around the user-selected locking grid point.

5. The method according to claim 4, characterized in that, The method further includes: When constructing the parameter boundaries of the objective function, in response to the locked grid points and locked values ​​selected by the user on the MAP graph, the upper and lower bounds of the parameter elements corresponding to the locked grid points are set within the tolerance range of the locked values, so as to apply an immovable hard constraint to the locked grid points.

6. The method according to claim 5, characterized in that, The method further includes: During the iterative optimization of the MAP grid node parameters, after each parameter update and when the final result is output, the parameter value corresponding to the user-selected locked grid point is forcibly restored to the locked value.

7. The method according to claim 1, characterized in that, The penalty items include smoothing penalty items; The smoothing penalty term is composed of the weighted sum of squared second-order differences in the X direction, the weighted sum of squared second-order differences in the Y direction, and a cross term, and is used to constrain the smoothness of the MAP surface.

8. The method according to claim 1, characterized in that, The penalty term includes a gradient constraint penalty term; The gradient constraint penalty term applies a unidirectional penalty to the first-order difference that exceeds the preset upper and lower limits of the gradient, which is used to constrain the rate of change of MAP in the coordinate axis direction.

9. The method according to claim 2, characterized in that, Based on the loss function value and the analytical gradient vector, an analytical gradient-based optimization algorithm with boundary constraints is used to iteratively optimize the MAP grid node parameters. In each iteration, boundary active set management is used to constrain parameter components that touch the boundary to the boundary values. Parameter updates are performed in the free parameter subspace to obtain the optimized MAP grid node parameters, including: Step a: In the current iteration step, determine the gradient direction based on the analyzed gradient vector; Step b: Based on the gradient direction, in each iteration, the parameter components that touch the boundary are constrained to the boundary value by the boundary active set management, and quasi-Newton optimization is performed in the free parameter subspace to determine the search direction and step size. Step c: Update the MAP grid node parameters according to the search direction and step size; Repeat steps a-c above until the convergence condition is met to obtain the optimized MAP grid node parameters.

10. The method according to claim 1, characterized in that, The method further includes: During the iterative optimization process, the first round starts from the original MAP value, performs a complete iterative optimization, and records the obtained optimal solution as the safety net benchmark. In subsequent rounds, starting from the current optimal solution, the MAP grid is downsampled to a low resolution, random noise is applied to the low-resolution grid, and then the original resolution is restored through interpolation upsampling, generating multiple perturbation-based candidate starting points; The optimal error is selected from the multiple candidate perturbation starting points by performing forward propagation evaluation, and the actual starting point is then re-executed from the actual starting point for a complete iterative optimization. If the results obtained from subsequent iterations are worse than the safety net benchmark and the degree of deterioration exceeds a first preset threshold, then the process reverts to the safety net benchmark and re-perturbs.

11. The method according to claim 10, characterized in that, The method further includes: If the error obtained from multiple consecutive optimization rounds does not decrease compared to the previous round, then increase the perturbation rate; If the error in this round decreases by more than the second preset threshold compared to the previous round, then the perturbation rate is reduced.

12. The method according to claim 1, characterized in that, Before constructing the objective function, the method further includes: Receive user-inputted custom fitting expression, parse the syntax structure of the expression, extract MAP lookup function calls, data column dependencies, and dependencies between Function nodes, and construct a computational dependency topology in the form of a directed acyclic graph; Pre-compute the non-dependent nodes in the computation-dependent topology, and retain the expression strings for the dependent nodes in the computation-dependent topology for subsequent dynamic evaluation; By replacing the discontinuous branch selection function in the expression with a continuously differentiable transition function, the computational graph of the traceable gradient propagation path is obtained.

13. The method according to claim 12, characterized in that, The method further includes: For each MAP lookup function in the computation graph, the interpolation weight and index of each data point are pre-calculated, and the interpolation weight and index are registered as a fixed buffer for lookup calls during the forward propagation process.

14. An engine MAP fitting and optimization device, characterized in that, include: A function construction module is used to construct an objective function, which includes an error term and at least one penalty term for constraining the characteristics of the MAP parameters; The computation module is used to perform one forward propagation and one backward propagation on the objective function based on a pre-built computation graph through an automatic differentiation framework, to obtain the loss function value of the objective function and the analytical gradient vector of the objective function with respect to all MAP grid node parameters; The iterative optimization module is used to iteratively optimize the MAP grid node parameters based on the loss function value and the analysis gradient vector using an analysis gradient-based optimization algorithm with boundary constraints. In each iteration, the parameter components that touch the boundary are constrained to the boundary value by using boundary active set management, and the parameter update is performed in the free parameter subspace to obtain the optimized MAP grid node parameters.

15. An electronic device, characterized in that, include: Memory, processor; The memory stores computer-executed instructions; The processor executes computer execution instructions stored in the memory to implement the method as described in any one of claims 1 to 13.

16. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer-executable instructions, which, when executed by a processor, are used to implement the method as described in any one of claims 1 to 13.

Images