Intelligent management method for large fixed-wing unmanned aerial vehicle flexible production line

Abstract

The present application relates to a kind of large fixed-wing unmanned aerial vehicle flexible production line's intelligent management method, including the digital twin scene of constructing large fixed-wing unmanned aerial vehicle flexible production line, digital twin scene includes multi-source heterogeneous data, multi-source heterogeneous data includes the digital model of each assembly station, logistics transport equipment, tooling and unmanned aerial vehicle product structure, the equipment operating state data of each assembly station, the position and operating parameter of logistics transport equipment, tooling usage state data and product assembly progress data are collected in real time by sensor network arranged in production site, it is associated with the corresponding digital model according to data type and collection position with binding;The present application is driven by digital twin, and production line state is evaluated and quantified in real time, combined with dynamic adaptability scheduling, self-adaptive rescheduling trigger and multi-objective assembly compensation optimization, the collaborative promotion of large fixed-wing unmanned aerial vehicle flexible mixed flow assembly efficiency and quality is realized.

Description

Technical Field

[0001] This invention relates to the field of production management, specifically to an intelligent management method for a flexible production line for large fixed-wing unmanned aerial vehicles. Background Technology

[0002] Large fixed-wing UAVs are characterized by their large structural dimensions, weak component rigidity, and high assembly precision requirements. Their production mode is shifting from mass production of a single model to mixed-flow production of multiple configurations and variable batches, placing extremely high demands on the flexibility and intelligent management capabilities of production lines. In existing technologies, flexible production lines face the following main technical challenges in mixed-flow assembly scenarios for multiple configuration UAVs: First, the dynamic evaluation methods for workstation status are limited, typically relying on equipment busy / idle conditions or single quality indicators for scheduling decisions. This fails to comprehensively reflect the impact of multiple dimensions such as tooling integrity, real-time load saturation, historical pass rates, and failure intervals on task-bearing capacity. This leads to a mismatch between task allocation and actual workstation capacity, resulting in both localized overload and idleness, hindering overall efficiency improvement. Second, there is a lack of quantitative analysis methods for the global impact of abnormal disturbances on the production line. When a delay occurs at a workstation, managers struggle to predict the cumulative effect of this delay propagating downstream along the logistics path and process sequence. Relying on experience to determine whether to intervene can easily lead to over- or under-response issues, affecting both production line stability and order delivery timeliness. Thirdly, during tooling changes or product configuration switching caused by rescheduling, the rapid adaptive compensation capability of the assembly baseline is insufficient. Traditional methods relying on repeated manual calibration are time-consuming and their accuracy is greatly affected by operator experience, making it difficult to meet the high-cycle, high-precision changeover requirements of flexible production lines. Fourthly, dynamic rescheduling decisions lack a multi-objective coordination mechanism that considers on-time delivery, load balancing, and adjustment stability. Existing methods are mostly single-objective optimizations or simple weighting, making it difficult to achieve a scientific balance among multiple conflicting objectives, resulting in rescheduling schemes failing to address all aspects in actual implementation. These problems collectively restrict the overall operating efficiency, anomaly response speed, and product assembly consistency of large fixed-wing UAV flexible production lines. To address these issues, an intelligent management method for large fixed-wing UAV flexible production lines is proposed. Summary of the Invention

[0003] To address the shortcomings of existing technologies, this invention provides an intelligent management method for a flexible production line of large fixed-wing UAVs, comprising the following steps: Step S1: Construct a digital twin scenario of a large fixed-wing UAV flexible production line. The digital twin scenario includes digital models of each assembly station, logistics and transportation equipment, tooling, and UAV product structure. Through a sensor network deployed on the production site, real-time collection of multi-source heterogeneous data such as equipment operation status data of each assembly station, location and operation parameters of logistics and transportation equipment, tooling usage status data, and product assembly progress data is carried out. The multi-source heterogeneous data is associated and bound with the corresponding digital model according to data type and collection location. According to the preset data mapping rules, the status parameters of physical entities are synchronized to the digital model in real time, driving the digital model to update the tooling in-situ status of each assembly station, the spatial pose of logistics and transportation equipment, the wear degree of tooling, and the assembly completion degree of UAV product structure, forming a dynamic production line status mirror that reflects the actual operating status of the current production line. Step S2: Based on the dynamic production line status mirror, and combined with the historical production efficiency and quality data of each workstation, calculate the dynamic adaptability index of each workstation for the current task type to be assigned. Step S3: Receive production orders containing drone products of different configurations, and generate an initial mixed-flow production line scheduling sequence that takes into account assembly priority and production line balance based on the process constraints of each configuration and the dynamic adaptability index of each workstation. Step S4: During the execution of the initial mixed-flow production line scheduling sequence, a directed graph of production line disturbance propagation is constructed with workstations as nodes and material flow and process sequence as directed edges. The cascading disturbance propagation intensity of each node is calculated in real time, thereby assessing the global blockage risk of the production line and generating a rescheduling trigger factor. Step S5: When the rescheduling triggering factor meets the preset triggering conditions, trigger dynamic rescheduling based on multi-objective optimization to generate a compensation scheme that includes the adjusted task allocation sequence and assembly baseline adaptive compensation parameters. Step S6: Perform simulation and verification of the compensation plan in the digital twin scenario. After successful verification, send it to the physical production line for execution and update the dynamic production line status image.

[0004] Furthermore, the construction of the digital twin scenario in step S1 specifically includes: Geometric model construction: Using the 3D CAD data from the production line design phase, 3D geometric models of each assembly station, logistics and transportation equipment, tooling, and drone product structure are established, and spatial assembly is carried out according to the actual production line layout to form a production line geometric scene. Physical property definition: Assign physical properties to each geometric model, including the kinematic parameters of logistics and transportation equipment, the travel range and repeatability of the positioners of each tooling, and the density and stiffness parameters of the materials of UAV components; Data association mapping: Establish an association mapping table between sensor data sources and digital models, bind the identifiers of each sensor to the attribute fields of the corresponding digital models, and define the acquisition frequency, data format conversion rules, and effective range verification rules for each sensor data. Communication interface configuration: A two-way data channel is established between the digital twin scenario and the various execution devices and measurement systems of the physical production line through the production line IoT platform, so that the digital twin scenario can receive positioning deviation data from online measurement feedback and issue compensation plans and scheduling instructions to the production line execution system.

[0005] Furthermore, in the calculation of the dynamic fitness index in step S2, it is necessary to normalize the indicators of each workstation. The specific normalization process is as follows: for each workstation, obtain the original value of the current tooling condition integrity of that workstation. Original value of load saturation And the original value of the historical first-piece pass rate extracted from historical production efficiency and quality data. and the original value of the number of times the historical average fault interval was handled. ; Among them, the original value of load saturation It is obtained by dividing the total standard working hours of the tasks currently awaiting processing at the workstation by the rated available working hours of the workstation within the statistical period. The original value of the current tooling condition integrity Original value of load saturation Historical first-piece pass rate (original value) and the original value of the number of times the historical average fault interval was handled. The normalized values ​​of tooling condition integrity were obtained by normalizing the values ​​using the range standardization method. Normalized value of load saturation Normalized value of first-piece pass rate and the normalized value of the number of fault interval handling times ; The range standardization process is as follows: For any original value X of an indicator, obtain the historical minimum value X of that indicator among similar workstations on the production line. min Compared to the historical maximum value X max The normalized value is obtained by subtracting the minimum value from the original value and then dividing by the difference between the maximum and minimum values. The original value X of any index is the original value of the current tooling condition integrity. Original value of load saturation Historical first-piece pass rate (original value) and the original value of the number of times the historical average fault interval was handled. Specifically: ; Among them, the normalized value of load saturation The process is reversed when calculating the dynamic fitness index, taking 1- This is to enable high adaptability to low load.

[0006] Furthermore, the dynamic fitness index is obtained by weighted fusion of normalized values. The specific calculation process is as follows: the normalized value M of the tooling condition integrity is... i Multiply by the first weighting coefficient α, and then calculate the inverted load saturation value 1-U. i Multiply by the second weighting coefficient β to normalize the first-piece yield rate to Q. i Multiply by the third weighting coefficient γ to normalize the number of fault interval processing values ​​F. i Multiply by the fourth weighting coefficient δ, and then add the four products together to obtain the dynamic fitness index of workstation i, which can be expressed by the formula: ; In the formula, α, β, γ and δ are all greater than 0 and their sum is 1; The weighting coefficients are determined by the analytic hierarchy process (AHP). Specifically, based on the relative importance of each indicator to the overall assembly capability of the workstation, a judgment matrix is ​​constructed using the 1-9 scale method. The eigenvector corresponding to the largest eigenvalue of the judgment matrix is ​​calculated and normalized to obtain the weighting coefficients. The weights are confirmed to be effective when the consistency ratio (CR) is less than 0.1.

[0007] Furthermore, the construction method of the directed graph for production line disturbance propagation in step S4 is as follows: each workstation is abstracted as a node, and the dependency relationship between workstations determined by the path of the automated guided vehicle and the process flow is abstracted as a directed edge, with the edge direction pointing from the immediate preceding workstation to the immediate following workstation; the edge weight w ij w represents the sensitivity of the completion delay of workstation i to the start time of its immediate successor workstation j. ij The calibration is performed by using the regression coefficient of the delay of workstation i and the waiting time of workstation j in historical production line operation data, or by measuring the offset of the start time of workstation j under different delay values ​​through digital twin simulation. The calculation of the cascaded disturbance transmission intensity includes: when any station i experiences an anomaly and generates an initial delay time. At that time, the delay propagates along the directed edge to the subsequent workstation j. The delay increment after propagation is calculated by multiplying the upstream delay by the edge weight and then subtracting the preset dynamic buffer time B of the downstream workstation. j And compare it with zero and take the larger value, specifically: ; In the formula, B j The dynamic buffer time for workstation j is determined by the slack time of the in-process task at that workstation. A rule was set that if multiple paths lead to the same workstation, the maximum value of the delay increment calculated for each path would be taken as the disturbance to that workstation. According to the above rules, the disturbance propagation proceeds downstream in the directed graph from the source workstation until it reaches the workstation with a zero out-degree in the directed graph of the production line disturbance propagation. This workstation is the end workstation of the production line. The delay increments on all propagation paths to the end workstation are accumulated to obtain the cascade disturbance propagation intensity C corresponding to workstation i. i ; When no abnormality occurs at workstation i, the cascade disturbance transmission intensity C i Set to zero; when multiple workstations experience anomalies simultaneously, calculate the cascade disturbance transmission intensity independently for each workstation with the abnormality as the source.

[0008] Furthermore, the rescheduling trigger factor is generated by obtaining the cascading disturbance propagation intensity C of all workstations. i、 Order urgency adjustment coefficient associated with each workstation Workstation importance weight For each workstation, its importance weight, urgency correction coefficient, and cascading disturbance propagation strength are multiplied together. The products of all workstations are then summed. The result is used as the rescheduling trigger factor. ; In the formula, N is the total number of workstations; Among them, the order urgency correction coefficient The calculation process is as follows: If workstation i currently has no orders in production, then the order urgency adjustment coefficient is... Take zero; If there are orders in production, the difference between the delivery deadline of the order and the current system time is taken as the remaining available time. The total time required for the unfinished assembly operations of the order at workstation i is taken as the required processing time. The ratio R of the remaining available time to the required processing time is calculated. i ; Obtain the preset maximum ratio threshold R in the system. max Subtract the ratio R of available time to required processing time from 1. i The maximum ratio threshold R max The ratio is used to obtain the order urgency correction factor. ,Right now ; Workstation Importance Weight The determination method is as follows: calculate the average work-in-process dwell time and process output rate of each station in the production line, identify the station with the longest average work-in-process dwell time and the lowest process output rate as the bottleneck station, and assign higher station importance weights to the bottleneck station and its immediate predecessor station. The value assigns a lower workstation importance weight to non-bottleneck workstations that have parallel substitution capabilities, and all weights are normalized to a sum of 1. The preset trigger condition is: when the rescheduling trigger factor... Greater than the adaptive threshold When this happens, dynamic rescheduling is initiated; The adaptive threshold is calculated as follows: a time window of length L is set, which slides forward with each system iteration, with a preset step size. Obtain all historical rescheduling trigger factor samples within the window and calculate the sample mean. with sample standard deviation The threshold is obtained by adding the product of the sensitivity coefficient k and the standard deviation to the mean. The sensitivity coefficient k is preset based on the production line's tolerance to disturbances, and its value ranges from 2 to 3.

[0009] Furthermore, the assembly datum adaptive compensation parameters in step S5 are generated in the following way: When rescheduling leads to tooling replacement or product configuration switching, the positioning deviation data of key feature points of the component to be assembled is obtained from the feedback of the online measurement system. The current working condition parameters and positioning deviation data are input into the pre-built compensation prediction model to obtain the initial compensation amount. In a digital twin scenario, an initial compensation amount is loaded, and a mechanism model incorporating assembly stress and geometric interference verification is used for verification. The root mean square value E of the residual at key points after compensation is then calculated. rms ; The calculation process for the root mean square value of the residual is as follows: obtain the deviation e between the compensated positions and the theoretical positions of all key feature points. j To calculate the square of each deviation, obtain the squared value, calculate the average of the squared values, and then take the square root. Specifically: ; In the formula, m is the number of key feature points; If the root mean square value of the residual at the key point after compensation is E rms If the value is less than the preset acceptance threshold, the initial compensation amount is used as the assembly datum adaptive compensation parameter; otherwise, based on the verification deviation feedback, the gradient descent algorithm is used to iteratively correct the compensation amount, using the root mean square value E of the residual at the key points after compensation. rms The loss function is defined by the compensation components of each positioner. To optimize the variables, the compensation amount is updated along the negative gradient direction. The update step size is preset based on the locator resolution and constrained by the physical travel limit of the locator. The iteration continues until the root mean square value E of the keypoint residuals after compensation is obtained. rms If the value is less than the acceptance threshold or the preset maximum number of iterations is reached, the final displacement compensation command value is output.

[0010] Furthermore, the compensation prediction model is a deviation compensation model based on the configuration adaptive influence coefficient matrix, and its construction and usage process includes: During the production line debugging phase, for each UAV product configuration and tooling combination, virtual assembly simulation is performed through a digital twin scenario. Unit displacement excitation with a preset amplitude is applied to each positioner in sequence. The preset amplitude is selected within the repeatability accuracy range of the positioner. The displacement response of the corresponding key feature points is recorded to form the excitation-response data pair of the combination. For each configuration and tooling combination, virtual simulations were repeated multiple times, and valid data recorded in the actual assembly history were superimposed. The data was then weighted and fused according to the weights of the virtual simulation data and the actual data. After removing samples whose deviations exceeded a preset anomaly threshold, the unit excitation vectors from each simulation were arranged column-wise into an excitation matrix. The response matrix R is formed by columnar arrangement of the displacement response vectors of the corresponding key points, and the influence coefficient matrix S under this combination is identified by weighted least squares method. c Its calculation method is to minimize the weighted sum of squared residuals, and the weight matrix used in the identification process is denoted as... Based on the preset measurement confidence level, specifically: ; Wherein, the subscript c represents the configuration-tooling combination index, and the superscript T represents the matrix transpose; During online forecasting, the corresponding influence coefficient matrix S is retrieved from the influence coefficient matrix library based on the product configuration and tooling used for the current task. c The key feature point positioning deviation vector obtained from online measurement As input, the compensation vector Construct a weighted objective function for the variables to be optimized. Where W is a positive definite weight matrix preset according to the assembly accuracy requirements of key feature points, and the superscript T indicates the matrix or vector transpose. Make the weighted objective function J with respect to The derivative is zero, and the optimal compensation vector obtained by solving is denoted as... Specifically: ; The optimal compensation vector This is the initial compensation amount.

[0011] Furthermore, the optimization objectives of the multi-objective optimization in step S5 include: minimizing the order delay penalty value, minimizing the variance of the working hours load of each workstation, and minimizing the total magnitude of the assembly baseline adaptive compensation parameters; The method for minimizing the order delay penalty value is as follows: For each order, subtract the delivery deadline from its estimated completion time; if the result is positive, take the difference; otherwise, take zero; then multiply by the unit time delay penalty weight corresponding to that order. Finally, sum the products of all orders as described above, expressed by the formula: ; In the formula, O represents the total number of orders, and C... o With D o These are the estimated completion time and delivery deadline for order o, respectively; The calculation method for the time load variance of each workstation is as follows: obtain the planned total processing time L for each workstation after rescheduling. i and the rated available working hours H of the workstation i Calculate the load rate of each workstation. Then calculate the load factor l of all workstations. i The variance; The total magnitude of the assembly datum adaptive compensation parameters is the final compensation parameter vector. The sum of the absolute values ​​of each component, i.e. ,in This is the compensation amount for the r-th positioner; The constraints for multi-objective optimization include at least the following: each workstation can only process one product at a time; the product process must follow the preset process sequence; the product start time must not be earlier than the material placement time; and the compensation amount of each positioner does not exceed the upper and lower limits of its physical stroke. The compensation scheme generation process adopts a hierarchical sequence method. Minimizing the order delay penalty P is given as the highest priority. The optimal value of this priority is solved, and all solutions whose target value deviation does not exceed the first preset tolerance are retained to form the current feasible solution set. Within the current feasible solution set, minimizing the variance of the working hours and load of each workstation is searched as the second highest priority. Solutions that satisfy this second-highest value are retained to form the updated feasible solution set. Finally, in the updated feasible solution set where the first two priorities are optimal, minimizing the total amplitude of the assembly benchmark adaptive compensation parameter is given as the final priority. The final scheduling sequence and compensation parameters are then obtained.

[0012] The present invention has the following advantages over the prior art: By constructing a digital twin scenario and collecting multi-source heterogeneous data from the production line in real time to form a dynamic state mirror, and based on this, by comprehensively considering multi-dimensional indicators such as tooling integrity, load saturation, first-piece pass rate, and number of fault interval handling times, the dynamic fitness index of the workstation is calculated using range standardization and analytic hierarchy process (AHP) weighted fusion. This accurately reflects the real-time undertaking capacity of each workstation for different configuration tasks, providing a precise quantitative basis for mixed-flow scheduling. Furthermore, by abstracting the material flow and process dependencies between workstations into a directed graph of disturbance propagation, and based on the edge weights calibrated by historical regression or simulation, and dynamic... The buffer time, calculating the cascading disturbance propagation intensity step by step along the directed edge, can accurately quantify the global impact of local anomalies on the end-of-line delivery capability. Combining the order urgency correction coefficient and the workstation importance weight determined by work-in-process dwell time and output rate, a rescheduling trigger factor is generated. An adaptive trigger threshold is constructed using the mean and standard deviation of historical samples within a sliding window, ensuring that the rescheduling timing precisely matches the actual risk level of the production line, avoiding both premature disruption to stable production and delayed severe congestion. When rescheduling triggers tooling changes or configuration switching, the... Using an influence coefficient matrix based on configuration-tooling combinations, a model is identified by virtual simulation excitation and weighted fusion of historical data for each configuration. This model rapidly converts the positioning deviation of key feature points measured online into initial compensation amounts. Furthermore, through iterative verification and gradient descent optimization using a digital twin mechanism model, the final compensation command value that minimizes the root mean square of the key point residuals is obtained. This significantly improves the assembly accuracy during rapid production changeover on flexible production lines while ensuring that assembly stress and geometric interference meet constraints. Simultaneously, with the optimization objectives of minimizing order delay penalties, minimizing workstation time load variance, and minimizing the total magnitude of assembly baseline compensation, a hierarchical sequence method is used to solve for the optimal solutions at each level under multiple hard constraints such as workstation exclusivity, process sequence, material placement, and positioner travel, while retaining the tolerance solution set. This approach balances order delivery timeliness, production line load balance, and adjustment operation stability. Finally, after verification through digital twin simulation, the system is deployed for execution, forming a complete closed-loop intelligent control system for flexible mixed-flow assembly of large fixed-wing UAVs. This effectively improves the overall operating efficiency of the production line, the speed of abnormal response, and the consistency of product assembly. Attached Figure Description

[0013] Figure 1 This is the overall flowchart of the present invention. Detailed Implementation

[0014] To further illustrate the technical means and effects of the present invention in achieving its intended purpose, the following detailed description of the specific implementation methods, structures, features, and effects of the present invention, in conjunction with the accompanying drawings and preferred embodiments, is provided below.

[0015] like Figure 1 As shown, an intelligent management method for a flexible production line of large fixed-wing UAVs includes the following steps: Step S1: Construct a digital twin scenario of a large fixed-wing UAV flexible production line. The digital twin scenario includes digital models of each assembly station, logistics and transportation equipment, tooling, and UAV product structure. Through a sensor network deployed on the production site, real-time collection of multi-source heterogeneous data such as equipment operation status data of each assembly station, location and operation parameters of logistics and transportation equipment, tooling usage status data, and product assembly progress data is carried out. The multi-source heterogeneous data is associated and bound with the corresponding digital model according to data type and collection location. According to the preset data mapping rules, the status parameters of physical entities are synchronized to the digital model in real time, driving the digital model to update the tooling in-situ status of each assembly station, the spatial pose of logistics and transportation equipment, the wear degree of tooling, and the assembly completion degree of UAV product structure, forming a dynamic production line status mirror that reflects the actual operating status of the current production line. Step S2: Based on the dynamic production line status mirror, and combined with the historical production efficiency and quality data of each workstation, calculate the dynamic adaptability index of each workstation for the current task type to be assigned. Step S3: Receive production orders containing drone products of different configurations, and generate an initial mixed-flow production line scheduling sequence that takes into account assembly priority and production line balance based on the process constraints of each configuration and the dynamic adaptability index of each workstation. Step S4: During the execution of the initial mixed-flow production line scheduling sequence, a directed graph of production line disturbance propagation is constructed with workstations as nodes and material flow and process sequence as directed edges. The cascading disturbance propagation intensity of each node is calculated in real time, thereby assessing the global blockage risk of the production line and generating a rescheduling trigger factor. Step S5: When the rescheduling triggering factor meets the preset triggering conditions, trigger dynamic rescheduling based on multi-objective optimization to generate a compensation scheme that includes the adjusted task allocation sequence and assembly baseline adaptive compensation parameters. Step S6: Perform simulation and verification of the compensation plan in the digital twin scenario. After successful verification, send it to the physical production line for execution and update the dynamic production line status image.

[0016] The construction of the digital twin scenario in step S1 specifically includes: Geometric model construction: Using the 3D CAD data from the production line design phase, 3D geometric models of each assembly station, logistics and transportation equipment, tooling, and drone product structure are established, and spatial assembly is carried out according to the actual production line layout to form a production line geometric scene. Physical property definition: Assign physical properties to each geometric model, including the kinematic parameters of logistics and transportation equipment, the travel range and repeatability of the positioners of each tooling, and the density and stiffness parameters of the materials of UAV components; Data association mapping: Establish an association mapping table between sensor data sources and digital models, bind the identifiers of each sensor to the attribute fields of the corresponding digital models, and define the acquisition frequency, data format conversion rules, and effective range verification rules for each sensor data. Communication interface configuration: A two-way data channel is established between the digital twin scenario and the various execution devices and measurement systems of the physical production line through the production line IoT platform, so that the digital twin scenario can receive positioning deviation data from online measurement feedback and issue compensation plans and scheduling instructions to the production line execution system; Through a standardized four-step process of geometric model construction, physical attribute definition, data association mapping, and communication interface configuration, the physical entity of a large fixed-wing UAV flexible production line is mapped to a digital twin scene with high fidelity. This ensures that the structural dimensions, motion constraints, accuracy indicators, and real-time status of each element of the production line correspond precisely in virtual space. In the geometric model construction, LOD lightweight technology is used to dynamically reduce the surface area of ​​the model outside the field of view, ensuring that the frame rate of large-scale scene rendering is no less than 30fps. In the multi-sensor data fusion, Kalman filtering is introduced to synchronously denoise the data from the laser tracker, force sensor, and temperature sensor, improving the reliability of the measurement data. At the same time, in the communication interface configuration, an XMLSchema message structure is defined, and a 10Hz timing synchronous acquisition loop and a linear extrapolation compensation mechanism triggered by packet loss of 3 consecutive frames are set to effectively avoid model jumps caused by data interruption. This provides a realistic, reliable, and highly robust virtual verification environment for subsequent dynamic fitness assessment, disturbance propagation simulation, and assembly compensation verification, significantly improving the rapid response capability and one-time assembly quality of the flexible production line in the face of multi-configuration mixed-flow assembly.

[0017] In the flexible production line for large fixed-wing drones, the docking station between the central wing and the outer wing is used to construct its digital twin scenario.

[0018] In the geometric model construction phase, a 3D CAD model of the production line is imported, including the docking platform, four CNC positioners, a laser tracker, and digital models of the central wing and outer wing products. Spatial assembly is completed according to the actual workshop coordinates. The theoretical positions of the ball heads of the four positioners in the global coordinate system are P1(12000,-800,3200), P2(12000,800,3200), P3(14000,-800,3200)), and P4(14000,800,3200), respectively, in millimeters. Simultaneously, LOD lightweighting is applied to the out-of-view equipment models, compressing the number of triangles to 20% of the original model to ensure real-time rendering performance. In the physical property definition phase, the travel of each axis of the positioners is assigned. Repeatability 0.03mm, maximum speed 5mm / s, and imparting density to the central wing aluminum alloy material. and elastic modulus This data is used for subsequent assembly mechanics simulation. In the data association and mapping stage, a mapping table is established between the laser tracker's measurement points and the corresponding target points in the digital model, and the original spherical coordinates of the tracker are... Convert to Cartesian coordinates (x, y, z) in real time using the following formula: ; Simultaneously, a coordinate validity range verification rule is set: if the converted coordinates deviate from the theoretical value... If it is found to be a gross error, it will be discarded.

[0019] For example, the coordinates of the target sphere at the process joint of the outer wing front beam were obtained in a certain data collection. The measured coordinates, after transformation and fusion of Kalman filtering and temperature compensation, are (3906.3, -887.4, 2099.8). The deviation from the theoretical value (3906, -887, 2100) is calculated as follows: ; If the deviation is less than the threshold, the data is valid, and the corresponding measurement point in the digital model is immediately updated to that coordinate. During the communication interface configuration phase, the laser tracker data stream is subscribed to via the OPCUA protocol through the production line IoT platform and pushed cyclically at a 10Hz time sequence. Each frame of the message encapsulates the data ID, timestamp, and 3D coordinates using a predefined XML format. If no valid message is received for three consecutive frames, linear extrapolation compensation is automatically triggered to estimate the current value based on the last valid value and the rate of change, ensuring data continuity. Simultaneously, the verified compensation scheme is sent to the CNC positioner PLC via the MQTT protocol, with the instruction delay controlled within 50ms. The digital twin scenario constructed in this way can synchronize the physical workstation status in real time with precision, providing a highly reliable decision-making basis for subsequent mixed-flow production line scheduling and assembly baseline adaptive compensation.

[0020] In step S2, during the calculation of the dynamic fitness index, it is necessary to normalize the indicators of each workstation. The specific normalization process is as follows: for each workstation, obtain the original value of the current tooling condition integrity of that workstation. Original value of load saturation And the original value of the historical first-piece pass rate extracted from historical production efficiency and quality data. and the original value of the number of times the historical average fault interval was handled. ; Among them, the original value of load saturation It is obtained by dividing the total standard working hours of the tasks currently awaiting processing at the workstation by the rated available working hours of the workstation within the statistical period. The original value of the current tooling condition integrity Original value of load saturation Historical first-piece pass rate (original value) and the original value of the number of times the historical average fault interval was handled. The normalized values ​​of tooling condition integrity were obtained by normalizing the values ​​using the range standardization method. Normalized value of load saturation Normalized value of first-piece pass rate and the normalized value of the number of fault interval handling times ; The range standardization process is as follows: For any original value X of an indicator, obtain the historical minimum value X of that indicator among similar workstations on the production line. min Compared to the historical maximum value X max The normalized value is obtained by subtracting the minimum value from the original value and then dividing by the difference between the maximum and minimum values. The original value X of any index is the original value of the current tooling condition integrity. Original value of load saturation Historical first-piece pass rate (original value) and the original value of the number of times the historical average fault interval was handled. Specifically: ; Among them, the normalized value of load saturation The process is reversed when calculating the dynamic fitness index, taking 1- This allows for high adaptability under low load. By clearly defining four original indicators reflecting the real-time task-taking capacity of a workstation—tool condition integrity, load saturation, first-piece pass rate, and number of fault interval handling times—and their acquisition methods, and employing a unified range standardization formula, the original values ​​with varying dimensions and ranges are uniformly mapped to the dimensionless interval [0,1] using historical minimum and maximum values ​​of similar workstations extracted from historical production efficiency and quality data. This eliminates scale differences between indicators and makes the performance of different workstations directly comparable under multi-dimensional conditions. Simultaneously, the normalized load saturation value is reversed, allowing workstations with lower current loads to obtain higher fitness scores. Thus, through the synergistic effect of normalization and reversal, the suitability of a workstation for undertaking new tasks is accurately quantified, providing standardized, unbiased, and physically meaningful basic data for subsequent weighted fusion and scheduling decisions.

[0021] Three wing docking stations in a flexible production line for large fixed-wing UAVs.

[0022] The current scheduling cycle requires selecting one of these three workstations to handle the upcoming wing docking task. The system obtains the raw values ​​of current indicators for each workstation from the dynamic production line status mirror and extracts the historical minimum and maximum values ​​of each indicator among similar workstations from historical production efficiency and quality data. The data is shown in the table below: Among them, the original value of load saturation Taking workstation 1 as an example, the acquisition method is as follows: the total standard working hours for the tasks currently awaiting processing at this workstation are 68 hours, and the rated available working hours within the statistical period are 80 hours. The same calculation applies to workstations 2 and 3.

[0023] The original values ​​of the four indicators were normalized using the range standardization formula.

[0024] Workstation 1: ; ; ; ; Workstation 2: ; ; ; ; Workstation 3: ; ; ; ; After normalization, all four indicators are mapped to the [0,1] interval and are dimensionless, allowing for direct comparison. Based on this, the normalized load saturation value is reversed: the load saturation of workstation 1 is as high as 0.85, after normalization... =0.833, inverse value This indicates that its available production capacity is extremely low; The load saturation of station 2 is only 0.4, after normalization. =0.333, inverse value 1- =0.667 indicates ample remaining capacity. Reversal processing ensures that low-load workstations have an advantage in subsequent weighted fusion, thus prioritizing tasks to workstations with sufficient idle capacity during scheduling, effectively avoiding localized overload.

[0025] The dynamic fitness index is obtained by weighted fusion of normalized values. The specific calculation process is as follows: the normalized value M of the tooling condition integrity is... i Multiply by the first weighting coefficient α, and then calculate the inverted load saturation value 1-U. i Multiply by the second weighting coefficient β to normalize the first-piece yield rate to Q. iMultiply by the third weighting coefficient γ to normalize the number of fault interval processing values ​​F. i Multiply by the fourth weighting coefficient δ, and then add the four products together to obtain the dynamic fitness index of workstation i, which can be expressed by the formula: ; In the formula, α, β, γ and δ are all greater than 0 and their sum is 1; The weight coefficients are determined by the analytic hierarchy process (AHP). Specifically, based on the relative importance of each indicator to the overall assembly capability of the workstation, a judgment matrix is ​​constructed using the 1-9 scale method. The eigenvector corresponding to the largest eigenvalue of the judgment matrix is ​​calculated and normalized to obtain the weight coefficients. The weights are confirmed to be effective when the consistency ratio (CR) is less than 0.1. By employing the analytic hierarchy process (AHP) to assess the relative importance of tooling integrity, load saturation, first-piece yield, and number of fault interval handling operations to the overall assembly capability of a workstation, a judgment matrix is ​​constructed using a 1-9 scale and its scientific weights are determined through consistency checks. The normalized and inverted values ​​of the aforementioned four indicators are then weighted and fused to ensure that the dynamic fitness index comprehensively reflects the multidimensional state of the workstation. Load saturation, with its highest weight, plays a dominant role, ensuring that low-load workstations receive higher scores. Meanwhile, tooling integrity and quality performance also influence the final decision, avoiding the bias of decisions based on a single indicator. This provides an objective and quantifiable basis for the optimal selection of task workstations in mixed-flow scheduling, significantly improving production line load balance and first-piece yield.

[0026] The index values ​​for the three wing docking stations, after normalization and reversal processing, are shown in the table below: First, the weight coefficients of each indicator are determined using the analytic hierarchy process (AHP).

[0027] The four indicators are tooling condition (M), load saturation (U), first piece pass rate (Q), and number of fault interval handling (F).

[0028] Based on practical experience in operating flexible production lines, process experts compared the relative importance of the four indicators pairwise, using a 1-9 scale (1 indicates equal importance, 3 indicates slightly important, 5 indicates significantly important, 7 indicates very important, 9 indicates extremely important, and the reciprocal indicates a negative comparison) to construct the following judgment matrix: ; The meaning of this matrix is: load saturation (U) is slightly more important than tooling condition (M) (scale 3), significantly more important than first-piece pass rate (Q) (scale 5), and significantly more important than failure interval handling number (F) (scale 5). Tooling condition integrity (M) is as important as first piece yield (Q) (scale 1), and slightly more important than failure interval handling (F) (scale 3). The first-piece yield (Q) is between equal and slightly important as the number of failure intervals handled (F) (scale 2).

[0029] This matrix aligns with the logic of prioritizing production line balance and then focusing on quality and equipment status in flexible mixed-flow manufacturing.

[0030] The weight vector is calculated using the square root method. First, the product of the elements in each row of the judgment matrix is ​​calculated: Product in the first row: ; Product in the second row: ; Product in row 3: ; Product in row 4: ; Take the fourth root of each product: ; ; ; ; The weight vector is obtained by normalizing the above square root: ; ; ; ; ; Verify the weight sum: The weights were fine-tuned to β=0.569 so that the sum was 1, and the final weights were α=0.194, β=0.569, γ=0.154, and δ=0.083.

[0031] Calculate the largest eigenvalue of the judgment matrix Perform a consistency check. First, calculate the weighted sum vector. : ; Calculate the ratio of each component to its corresponding weight: ; ; ; ; The approximate value of the largest eigenvalue is: ; Consistency Indicators .

[0032] When n=4, the average random consistency index RI=0.9, and the consistency ratio .

[0033] Since CR = 0.030 < 0.1, the judgment matrix has satisfactory consistency, and the obtained weights are valid. According to the weighted fusion formula Calculate the dynamic adaptability index of each workstation one by one.

[0034] Workstation 1: ; Workstation 2: ; Workstation 3: ; The calculation results show that the dynamic adaptability index of station 2 is the highest (0.694). Although its tooling condition is slightly lower, it has the strongest overall acceptance capability due to its extremely low load saturation (contribution of up to 0.3795 after reversal) and excellent first-piece pass rate. Workstation 3 is in the middle; workstation 1 received the lowest score because its load was nearing saturation and its quality performance was mediocre. Based on this, the scheduling system prioritized assigning the next wing docking task to workstation 2, thereby effectively utilizing idle capacity, delaying the backlog of tasks at high-load workstations, and ensuring quality reliability, thus achieving scientific and efficient task allocation on the flexible production line.

[0035] The construction method of the directed graph for production line disturbance propagation in step S4 is as follows: each workstation is abstracted as a node, and the dependency relationship between workstations determined by the path of the automated guided vehicle and the process flow is abstracted as a directed edge, with the edge direction pointing from the immediate predecessor workstation to the immediate successor workstation; the edge weight w ij w represents the sensitivity of the completion delay of workstation i to the start time of its immediate successor workstation j. ij The calibration is performed by using the regression coefficient of the delay of workstation i and the waiting time of workstation j in historical production line operation data, or by measuring the offset of the start time of workstation j under different delay values ​​through digital twin simulation. The calculation of the cascaded disturbance transmission intensity includes: when any station i experiences an anomaly and generates an initial delay time. At that time, the delay propagates along the directed edge to the subsequent workstation j. The delay increment after propagation is calculated by multiplying the upstream delay by the edge weight and then subtracting the preset dynamic buffer time B of the downstream workstation. j And compare it with zero and take the larger value, specifically: ; In the formula, B j The dynamic buffer time for workstation j is determined by the slack time of the in-process task at that workstation. A rule was set that if multiple paths lead to the same workstation, the maximum value of the delay increment calculated for each path would be taken as the disturbance to that workstation. According to the above rules, the disturbance propagation proceeds downstream in the directed graph from the source workstation until it reaches the workstation with a zero out-degree in the directed graph of the production line disturbance propagation. This workstation is the end workstation of the production line. The delay increments on all propagation paths to the end workstation are accumulated to obtain the cascade disturbance propagation intensity C corresponding to workstation i. i ; When no abnormality occurs at workstation i, the cascade disturbance transmission intensity C i Set to zero; when multiple workstations experience anomalies simultaneously, calculate the cascaded disturbance transmission intensity independently for each workstation with the anomaly as the source. By abstracting the dependencies between workstations determined by the paths of automated guided vehicles and the flow of processes into directed edges, and using historical operational data regression or digital twin simulation to calibrate edge weights to quantify the sensitivity of delay propagation, and combining the dynamic buffer time determined by the relaxation time of each workstation based on the in-process task, the cascading disturbance propagation intensity is calculated in a step-by-step propagation and maximum value manner. This can accurately characterize the cumulative impact of local anomalies along the production line topology to the final delivery node, thereby accurately distinguishing the contribution of anomalies at different workstations to the global blockage risk. This provides precise quantitative input for the generation of rescheduling trigger factors, avoiding a one-size-fits-all over-response to all anomalies or an insufficient response to critical bottleneck anomalies, and effectively ensuring the stability and on-time delivery capability of flexible production lines when facing disturbances in multi-configuration mixed-flow assembly.

[0036] A section of the wing assembly serial production line in a large fixed-wing UAV flexible production line contains four workstations, with the process sequence as follows: Workstation 1 (wing rib assembly) → Workstation 2 (skin riveting) → Workstation 3 (system installation) → Workstation 4 (final inspection and delivery).

[0037] Material transfer between workstations is completed by automated guided vehicles, and the process flow dependencies are determined. Based on regression analysis of historical operating data and digital twin simulation verification, the weights of each edge are calibrated as follows: w 12 =0.9 (Sensitivity of the effect of station 1 delay on station 2), w 23 =0.8, w 34 =0.7.

[0038] Each workstation presets a dynamic buffer time based on the relaxation time of the currently in-process task: .

[0039] At the current moment, an anomaly has occurred at workstation 1 due to a tooling malfunction, and an initial delay is expected. .

[0040] According to the calculation method for the conduction strength of cascaded disturbances, the propagation process is as follows: Workstation 1 → Workstation 2: ; Workstation 2 → Workstation 3: ; Workstation 3 → Workstation 4: ; The end station of the production line is station 4, with an out-degree of zero. Accumulating the delay increment along the unique propagation path to the end station yields the cascaded disturbance propagation intensity corresponding to station 1. (because ).

[0041] If the calculated end-point delay increment is non-zero, then this value should be accumulated. Here, the end-point delay is absorbed by the buffer, indicating that the initial 20-minute delay at workstation 1 does not substantially delay the final delivery under the existing buffer, and the cascading risk is low. If an anomaly occurs at workstation 2 simultaneously, the initial delay... Then, calculate its impact on the end independently: first propagate to workstation 3, Then proceed to workstation 4. .

[0042] If an abnormality occurs at workstation 3 ,but Therefore .

[0043] It is evident that the anomaly at workstation 3 directly penetrates the buffer, causing significant delay at the end point. The cascading disturbance propagation intensity is the strongest, and this should be given special attention. This method is a rescheduling trigger factor. It provides accurate workstation-level risk contribution values, avoiding potential misjudgments that might arise from judging global risks based solely on the timing of single-point anomalies.

[0044] The rescheduling trigger factor is generated by obtaining the cascaded disturbance propagation intensity C of all workstations. i、 Order urgency adjustment coefficient associated with each workstation Workstation importance weight For each workstation, its importance weight, urgency correction coefficient, and cascading disturbance propagation strength are multiplied together. The products of all workstations are then summed. The result is used as the rescheduling trigger factor. ; In the formula, N is the total number of workstations; Among them, the order urgency correction coefficient The calculation process is as follows: If workstation i currently has no orders in production, then the order urgency adjustment coefficient is... Take zero; If there are orders in production, the difference between the delivery deadline of the order and the current system time is taken as the remaining available time. The total time required for the unfinished assembly operations of the order at workstation i is taken as the required processing time. The ratio R of the remaining available time to the required processing time is calculated. i ; Obtain the preset maximum ratio threshold R in the system. max Subtract the ratio R of available time to required processing time from 1. i The maximum ratio threshold R max The ratio is used to obtain the order urgency correction factor. ,Right now ; Workstation Importance Weight The determination method is as follows: calculate the average work-in-process dwell time and process output rate of each station in the production line, identify the station with the longest average work-in-process dwell time and the lowest process output rate as the bottleneck station, and assign higher station importance weights to the bottleneck station and its immediate predecessor station. The value assigns a lower workstation importance weight to non-bottleneck workstations that have parallel substitution capabilities, and all weights are normalized to a sum of 1. The preset trigger condition is: when the rescheduling trigger factor... Greater than the adaptive threshold When this happens, dynamic rescheduling is initiated; The adaptive threshold is calculated as follows: a time window of length L is set, which slides forward with each system iteration, with a preset step size. Obtain all historical rescheduling trigger factor samples within the window and calculate the sample mean. with sample standard deviation The threshold is obtained by adding the product of the sensitivity coefficient k and the standard deviation to the mean. The sensitivity coefficient k is preset based on the production line's tolerance to disturbances, and its value ranges from 2 to 3.

[0045] By multiplying and summing the cascading disturbance propagation intensity of each workstation, the order urgency correction coefficient, and the workstation importance weight, a quantitative rescheduling trigger factor is generated that comprehensively reflects the scope of disturbance impact, delivery urgency, and production line bottleneck location. The order urgency correction coefficient accurately characterizes the urgency of work-in-process orders by normalizing the ratio of remaining available time to required processing time. The workstation importance weight automatically identifies bottleneck workstations and their immediate predecessor workstations based on average work-in-process dwell time and process output rate, assigning differentiated importance to them. This allows the trigger factor to focus on the disturbance sensitivity of high-urgency orders and high-bottleneck workstations. At the same time, an adaptive threshold is constructed using the mean and standard deviation of historical trigger factors within a sliding time window, and the trigger sensitivity is adjusted by a sensitivity coefficient. This ensures that dynamic rescheduling is only triggered when disturbances may cause order delays or severe production line blockages. This avoids system oscillations caused by frequent scheduling and prevents delivery defaults caused by response lags, significantly improving the accuracy and stability of scheduling decisions for flexible mixed-flow production lines in disturbance environments.

[0046] The wing assembly serial production line in the flexible production line for large fixed-wing UAVs contains four workstations.

[0047] At the current moment, the production line status is as follows: both station 1 and station 2 have experienced anomalies, and the cascade disturbance transmission intensity of each station is calculated independently as C1=0, C2=0, C3=9.5, and C4=0 (station 4 is the end station and does not generate transmission intensity itself).

[0048] The urgency information for orders in production at each workstation is as follows: Workstation 1 currently has in-process order A. The difference between the delivery deadline and the current system time is the remaining available time. The total time required for the assembly process that has not yet been completed is Calculate R1 = 480 / 240 = 2.

[0049] Workstation 2 has order B, remaining available time Work hours required Calculate R² = 120 / 300 = 0.4.

[0050] Workstation 3 has order C, remaining available time Work hours required , .

[0051] Workstation 4 has no work-in-process orders (final inspection and delivery station, no urgency level), therefore The system's preset maximum ratio threshold R max =5.

[0052] Calculate the order urgency correction factor for each workstation. : ; ; ; (No work-in-process orders).

[0053] Workstation Importance Weight The bottleneck station is determined based on the average work-in-process dwell time and process output rate of each station. According to historical data statistics, station 3 (system installation) has the longest average work-in-process dwell time (up to 45 minutes) and the lowest process output rate (1.33 pieces / hour), and is therefore identified as the bottleneck station; its immediate preceding station is station 2 (skin riveting).

[0054] Based on this, the following weights are assigned: Bottleneck workstation 3 Immediately preceding station 2 Workstation 1 Workstation 4 is not a bottleneck and has no immediate preceding effects. All weights sum to 1.

[0055] Using formula Calculate the rescheduling trigger factor: ; Rescheduling trigger factor .

[0056] Now calculate the adaptive threshold. Set the time window length L = 30 minutes, and the sliding step size... The current window contains the 6 most recent historical trigger factor samples (units are the same as Γ): {1.2, 0.8, 1.5, 2.1, 1.9, 2.4}. Calculate the sample mean: ; Sample standard deviation: ; With a sensitivity coefficient k = 2.5 (within the range of 2 to 3), calculate the adaptive threshold: ; The comparison shows that the current triggering factor The triggering conditions have not been met, and the system will not initiate dynamic rescheduling for the time being. Although the cascading disturbance propagation intensity of workstation 3 is relatively high (9.5 minutes), and its order urgency coefficient and workstation importance weight are both large, considering the overall disturbance and historical fluctuations, the current comprehensive risk has not significantly exceeded the recent normal fluctuation range, indicating that the production line still has the ability to absorb disturbances through its own buffer, and no immediate intervention is required. If the delay of workstation 3 further intensifies, or the order urgency increases, If the threshold is exceeded, a rescheduling will be triggered, thus avoiding unnecessary scheduling adjustments that could impact production line stability.

[0057] In step S5, the assembly datum adaptive compensation parameters are generated in the following way: When rescheduling leads to tooling replacement or product configuration switching, the positioning deviation data of key feature points of the component to be assembled is obtained from the feedback of the online measurement system. The current working condition parameters and positioning deviation data are input into the pre-built compensation prediction model to obtain the initial compensation amount. In a digital twin scenario, an initial compensation amount is loaded, and a mechanism model incorporating assembly stress and geometric interference verification is used for verification. The root mean square value E of the residual at key points after compensation is then calculated. rms ; The calculation process for the root mean square value of the residual is as follows: obtain the deviation e between the compensated positions and the theoretical positions of all key feature points. j To calculate the square of each deviation, obtain the squared value, calculate the average of the squared values, and then take the square root. Specifically: ; In the formula, m is the number of key feature points; If the root mean square value of the residual at the key point after compensation is E rms If the value is less than the preset acceptance threshold, the initial compensation amount is used as the assembly datum adaptive compensation parameter; otherwise, based on the verification deviation feedback, the gradient descent algorithm is used to iteratively correct the compensation amount, using the root mean square value E of the residual at the key points after compensation. rms The loss function is defined by the compensation components of each positioner. To optimize the variables, the compensation amount is updated along the negative gradient direction. The update step size is preset based on the locator resolution and constrained by the physical travel limit of the locator. The iteration continues until the root mean square value E of the keypoint residuals after compensation is obtained. rms If the value is less than the acceptance threshold or the preset maximum number of iterations is reached, the final displacement compensation command value is output. When rescheduling triggers tooling changes or product configuration switching, the positioning deviation of key feature points is obtained through an online measurement system, and an initial compensation amount is generated using a compensation prediction model to quickly respond to deviation changes. The initial compensation amount is loaded in the digital twin scenario and verified with high fidelity using a mechanism model that includes assembly stress and geometric interference verification. The root mean square value of the residuals of the key points after compensation is used as a quantitative evaluation index of assembly accuracy to objectively reflect the current compensation effect. When the accuracy does not meet the requirements, a gradient descent algorithm is used with the root mean square value of the residuals as the loss function to iteratively correct the compensation amount of each positioner along the negative gradient direction. The update step size is preset according to the positioner resolution and strictly limited by the physical travel, ensuring that each correction is within the achievable accuracy and travel range of the equipment. Thus, under the premise of ensuring assembly stress and no interference constraints, the assembly deviation is gradually converged to within the acceptance threshold. Finally, a displacement compensation command value that can be directly issued and executed is output, which effectively improves the first-pass yield and adjustment efficiency of the flexible production line during rapid production changeover and avoids assembly quality fluctuations caused by tooling or configuration changes.

[0058] In the wing docking scenario, after rescheduling, it is necessary to switch from type A UAV to type B UAV, and at the same time change the docking fixtures for the middle and outer wings. An online laser tracker measures the positioning deviation vectors of two key feature points on the middle and outer wing components. .

[0059] The compensation prediction model quickly solves the influence coefficient matrix of this configuration and tooling combination, and outputs the initial compensation vector. This corresponds to the displacement adjustment command of two CNC positioners.

[0060] The initial compensation amount is loaded into the digital twin scenario, and virtual assembly is performed using a finite element mechanism model that includes assembly stress and geometric interference verification to obtain the residual deviation after compensation at two key points. Calculate the root mean square value of the residuals: ; Preset acceptance threshold ,because The initial compensation scheme does not meet the accuracy requirements, so gradient descent iterative correction is initiated.

[0061] by Let the loss function be the loss function at the current compensation point. At this point, by applying independent perturbations with a step size of 0.001 mm to the digital twin mechanism model, the Jacobian matrix is ​​numerically estimated: ; Calculate the gradient vector of the loss function : ; ; The preset update step size η = 0.005 mm (matching the positioner resolution) is used to update the compensation amount along the negative gradient direction. ; After the update, the compensation amounts for each locator are all within the range of... Within the travel range, the limits were not exceeded.

[0062] The corrected compensation amount is then reloaded into the digital twin mechanism model for verification, and this process is repeated iteratively until... Down to Once the acceptance criteria are met, the final displacement compensation command value is output for the physical locator to execute.

[0063] The entire process organically closes the loop between online measurement, rapid prediction, high-fidelity verification, and constrained iterative optimization, significantly improving the adaptive compensation accuracy and engineering feasibility of assembly references during configuration switching.

[0064] The compensation prediction model is a deviation compensation model based on the configuration adaptive influence coefficient matrix. Its construction and usage process includes: During the production line debugging phase, for each UAV product configuration and tooling combination, virtual assembly simulation is performed through a digital twin scenario. Unit displacement excitation with a preset amplitude is applied to each positioner in sequence. The preset amplitude is selected within the repeatability accuracy range of the positioner. The displacement response of the corresponding key feature points is recorded to form the excitation-response data pair of the combination. For each configuration and tooling combination, virtual simulations were repeated multiple times, and valid data recorded in the actual assembly history were superimposed. The data was then weighted and fused according to the weights of the virtual simulation data and the actual data. After removing samples whose deviations exceeded a preset anomaly threshold, the unit excitation vectors from each simulation were arranged column-wise into an excitation matrix. The response matrix R is formed by columnar arrangement of the displacement response vectors of the corresponding key points, and the influence coefficient matrix S under this combination is identified by weighted least squares method. c Its calculation method is to minimize the weighted sum of squared residuals, and the weight matrix used in the identification process is denoted as... Based on the preset measurement confidence level, specifically: ; Wherein, the subscript c represents the configuration-tooling combination index, and the superscript T represents the matrix transpose; During online forecasting, the corresponding influence coefficient matrix S is retrieved from the influence coefficient matrix library based on the product configuration and tooling used for the current task. c The key feature point positioning deviation vector obtained from online measurement As input, the compensation vector Construct a weighted objective function for the variables to be optimized. Where W is a positive definite weight matrix preset according to the assembly accuracy requirements of key feature points, and the superscript T indicates the matrix or vector transpose. Make the weighted objective function J with respect to The derivative is zero, and the optimal compensation vector obtained by solving is denoted as... Specifically:

[0065] The optimal compensation vector This is the initial compensation amount; For different UAV product configurations and tooling combinations, small-amplitude unit displacement excitations are sequentially applied to each locator in a digital twin scenario during the offline phase. Multiple virtual simulations and weighted historical assembly data are integrated, and after removing outliers, a weighted least squares method is used to identify the influence coefficient matrix specific to the configuration-tooling combination. The weighted fusion of virtual and real data effectively overcomes the shortcomings of pure simulation model bias and actual data noise, enabling the matrix to accurately reflect the linear influence of locator displacement on key feature points under this combination. During online prediction, the corresponding matrix is ​​quickly retrieved based on the current configuration and tooling, and the measurement deviation is mapped to the locator compensation amount in one go using a weighted pseudo-inverse method. The solution process considers the weight differences in the accuracy requirements of each key point, automatically coordinating the trade-offs of deviations at each point while satisfying least squares optimality. This provides a high-confidence initial compensation amount for subsequent digital twin mechanism verification, significantly improving the response speed and solution quality of assembly benchmark adaptive compensation during configuration switching and tooling replacement, effectively ensuring the positioning accuracy and adjustment efficiency of flexible production lines in multi-configuration mixed-flow assembly.

[0066] The outer wing docking station of the large fixed-wing UAV flexible production line needs to be adapted to the B-type UAV and its special tooling combination, which is denoted as configuration-tooling combination index c=1.

[0067] In the offline phase, unit displacement excitations are sequentially applied to two CNC positioners using digital twin virtual assembly simulation, with amplitudes of 0.1 mm for each (within the positioner's repeatability accuracy). (Within the scope). A total of two virtual simulation stimuli and one actual assembly data acquisition were performed. The data is as follows: Virtual incentive 1: Key point response ; Virtual Incentive 2: ,response ; Actual historical incentives: ,response .

[0068] Preset measurement confidence weight matrix (The virtual simulation data has a weight of 1, and the actual data has a weight of 0.5). After removing outlier samples, the excitation vectors and response vectors are arranged into an excitation matrix column by column. and response matrix Using the formula Identify the influence coefficient matrix.

[0069] Calculated .

[0070] During online prediction, before the docking of the outer wings of the Type B UAV, the positioning deviation vector of key feature points is obtained through online measurement. .

[0071] Preset positive definite weight matrix W=I (with the same accuracy requirement for each point), according to the formula Solve for the initial compensation amount: ; ; ; ; Will As the initial compensation amount, it is provided to the digital twin mechanism model for verification and iterative correction. After gradient descent iterative optimization, it quickly converges to the final compensation parameter that meets the assembly accuracy requirements. This realizes a complete closed loop from deviation measurement and rapid prediction of the influence coefficient matrix to high-fidelity verification and correction, effectively ensuring the adaptive compensation accuracy and engineering practicality of the assembly datum during configuration switching.

[0072] The optimization objectives of the multi-objective optimization in step S5 include: minimizing the order delay penalty value, minimizing the variance of the working hours load of each workstation, and minimizing the total magnitude of the assembly baseline adaptive compensation parameters; The method for minimizing the order delay penalty value is as follows: For each order, subtract the delivery deadline from its estimated completion time; if the result is positive, take the difference; otherwise, take zero; then multiply by the unit time delay penalty weight corresponding to that order. Finally, sum the products of all orders as described above, expressed by the formula: ; In the formula, O represents the total number of orders, and C... o With D o These are the estimated completion time and delivery deadline for order o, respectively; The calculation method for the time load variance of each workstation is as follows: obtain the planned total processing time L for each workstation after rescheduling. i and the rated available working hours H of the workstation i Calculate the load rate of each workstation. Then calculate the load factor l of all workstations. iThe variance; The total magnitude of the assembly datum adaptive compensation parameters is the final compensation parameter vector. The sum of the absolute values ​​of each component, i.e. ,in This is the compensation amount for the r-th positioner; The constraints for multi-objective optimization include at least the following: each workstation can only process one product at a time; the product process must follow the preset process sequence; the product start time must not be earlier than the material placement time; and the compensation amount of each positioner does not exceed the upper and lower limits of its physical stroke. The compensation scheme generation process adopts a hierarchical sequence method. Minimizing the order delay penalty value P is taken as the highest priority. The optimal value of this priority is solved and all solutions whose target value deviation does not exceed the first preset tolerance are retained to form the current feasible solution set. Within the current feasible solution set, minimizing the variance of the working hours and load of each workstation is taken as the second highest priority for searching. Solutions that satisfy this second-highest value are retained to form the updated feasible solution set. Finally, in the updated feasible solution set where the first two priorities are optimal, minimizing the total amplitude of the assembly benchmark adaptive compensation parameter is taken as the final priority to obtain the final scheduling sequence and compensation parameters. In dynamic rescheduling, three optimization objectives are considered simultaneously: minimizing the order delay penalty, the variance of the workload of each workstation, and the total magnitude of the assembly baseline adaptive compensation parameters. A hierarchical sequence method is used to solve these objectives step by step according to priority. First, the highest priority objective of on-time order delivery is optimized while retaining multiple feasible solutions within the tolerance range. Second, the scheduling scheme that best balances the workload of each workstation is selected from the candidate solutions. Finally, under the premise that the first two objectives have been satisfied, the total magnitude of the positioner compensation action is further minimized. This approach balances on-time delivery, production line balance, and the smoothness of adjustment operations while avoiding the defects of multi-objective weighted methods, such as the difficulty in determining weights and the mutual cancellation between objectives. It can generate an executable scheduling and compensation scheme that satisfies delivery commitments, maintains stable production line operation, and minimizes the impact on equipment when flexible mixed-flow assembly faces disturbances, significantly improving the overall quality and engineering practicality of rescheduling decisions.

[0073] In the wing docking scenario, a sudden anomaly at workstation 2 has triggered a rescheduling. The production line needs to reassign two pending orders (Order A and Order B) across three workstations (Workstation 1, Workstation 2, and Workstation 3). Simultaneously, tooling changes require determining the compensation amounts for two positioners. The known information is as follows: Order A: Delivery Deadline Delay penalty weight per unit time ; Order B: Delivery Deadline Delay penalty weight per unit time .

[0074] The rated available working hours for each workstation are The current time is t=0.

[0075] Based on process constraints and workstation adaptability, there are only two feasible scheduling options (which already satisfy hard constraints such as workstation exclusivity, process sequence, and material availability). The estimated completion time and total planned processing time for each workstation, as well as the compensation parameters, are as follows: Option 1: Assign order A to workstation 3 and order B to workstation 1.

[0076] Estimated completion time: .

[0077] Total planned processing time for each workstation: Workstation 1 (Processing Order B), Workstation 2 (Abnormal shutdown), workstation 3 (Processing order A).

[0078] The assembly datum compensation parameter vector corresponding to this scheme is: .

[0079] Option 2: Assign order A to workstation 1 and order B to workstation 3.

[0080] Estimated completion time: .

[0081] Total planned processing time for each workstation: Workstation 1 Workstation 2 Workstation 3 .

[0082] Compensation parameter vector .

[0083] First, calculate the highest priority objective, the order delay penalty value P: ; Option 1: ; Option 2: ; Both schemes have an optimal order delay penalty value of 0, which is within the tolerance range (the default tolerance is 0), and both are retained to proceed to the next priority level.

[0084] The second priority objective is the variance of work hours at each workstation. First, calculate the workstation load rate for each plan. .

[0085] Option 1 load rate: .

[0086] mean .

[0087] variance: ; Option 2 load rate: .

[0088] mean .

[0089] variance: ; Obviously Option 1 has a more balanced load, so it is better in the second priority category. Option 1 is retained and placed in the last priority category.

[0090] Third priority objective – Total magnitude of assembly datum adaptive compensation parameters .

[0091] The total compensation amount for Option 1 is: ; (Scheme 2 has been eliminated and does not require calculation.) Ultimately, the hierarchical sequence method outputs Scheme 1 as the final scheduling sequence, with a total compensation parameter amplitude of 0.8mm, small positioner movement, and smooth execution. The entire process ensures zero order delays, the most balanced load, and the most economical compensation actions, achieving multi-objective comprehensive optimization of rescheduling decisions under disturbances.

[0092] The above are merely preferred embodiments of the present invention and are not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. An intelligent management method for a flexible production line of large fixed-wing unmanned aerial vehicles, characterized in that, Includes the following steps: Step S1: Construct a digital twin scenario of a large fixed-wing UAV flexible production line. The digital twin scenario includes multi-source heterogeneous data, which includes digital models of each assembly station, logistics and transportation equipment, tooling, and UAV product structure. Through a sensor network deployed on the production site, the operating status data of each assembly station, the location and operating parameters of the logistics and transportation equipment, the usage status data of the tooling, and the product assembly progress data are collected in real time. The multi-source heterogeneous data are associated and bound with the corresponding digital models according to the data type and collection location. According to the preset data mapping rules, the status parameters of the physical entities are synchronized to the digital model in real time for processing to generate a dynamic production line status mirror of the current actual operating status of the production line. Step S2: Based on the dynamic production line status mirror, and combined with the historical production efficiency and quality data of each workstation, calculate the dynamic adaptability index of each workstation for the current task type to be assigned. Step S3: Receive production orders containing drone products of different configurations, and generate an initial mixed-flow production line scheduling sequence that takes into account assembly priority and production line balance based on the process constraints of each configuration and the dynamic adaptability index of each workstation. Step S4: During the execution of the initial mixed-flow production line scheduling sequence, a directed graph of production line disturbance propagation is constructed with workstations as nodes and material flow and process sequence as directed edges. The cascading disturbance propagation intensity of each node is calculated in real time, thereby assessing the global blockage risk of the production line and generating a rescheduling trigger factor. Step S5: When the rescheduling triggering factor meets the preset triggering conditions, trigger dynamic rescheduling based on multi-objective optimization to generate a compensation scheme that includes the adjusted task allocation sequence and assembly baseline adaptive compensation parameters. Step S6: Perform simulation and verification of the compensation plan in the digital twin scenario. After successful verification, send it to the physical production line for execution and update the dynamic production line status image.

2. The intelligent management method for a flexible production line of a large fixed-wing UAV according to claim 1, characterized in that: In the process of calculating the dynamic fitness index in step S2, it is necessary to normalize the indicators of each workstation. The specific normalization process is as follows: for each workstation, obtain the original value of the current tooling status integrity, the original value of the load saturation, the original value of the historical first piece pass rate extracted from historical production efficiency data and quality data, and the original value of the historical average fault interval processing times. The original value of load saturation is obtained by dividing the total standard working hours of the current tasks to be processed at the workstation by the rated available working hours of the workstation within the statistical period. The original values ​​of current tooling condition integrity, load saturation, historical first-piece pass rate, and historical average fault interval processing times are normalized using the range standardization method to obtain normalized values ​​of tooling condition integrity, load saturation, first-piece pass rate, and fault interval processing times. The range standardization process is as follows: For any original value of an indicator, obtain the historical minimum and historical maximum values ​​of the indicator in the same type of workstation on the production line. Subtract the minimum value from the original value and divide by the difference between the maximum and minimum values ​​to obtain the normalized value. The original values ​​of any indicator include the original value of the current tooling condition integrity, the original value of the load saturation, the original value of the historical first piece pass rate, and the original value of the historical average fault interval processing times. The load saturation normalization value is reversed when the dynamic fitness index is calculated in the subsequent calculation. The load saturation normalization value is taken as 1 and subtracted from the load saturation normalization value so that low load corresponds to high fitness.

3. The intelligent management method for a flexible production line of a large fixed-wing UAV according to claim 2, characterized in that: The dynamic fitness index is obtained by weighted fusion of normalized values. The specific calculation process is as follows: multiply the normalized value of tooling condition integrity by the first weight coefficient, multiply the value of load saturation after reversal by the second weight coefficient, multiply the normalized value of first piece pass rate by the third weight coefficient, multiply the normalized value of fault interval processing times by the fourth weight coefficient, and then add the four products to obtain the dynamic fitness index of workstation i. The weight coefficients are determined by the analytic hierarchy process (AHP). Specifically, based on the relative importance of each indicator to the overall assembly capability of the workstation, a judgment matrix is ​​constructed using the 1-9 scale method. The eigenvector corresponding to the largest eigenvalue of the judgment matrix is ​​calculated and normalized to obtain the weight coefficients. The weights are confirmed to be valid when the consistency ratio is less than the preset value.

4. The intelligent management method for a flexible production line of a large fixed-wing UAV according to claim 3, characterized in that: The construction method of the directed graph of production line disturbance transmission in step S4 is as follows: each workstation is abstracted as a node, and the dependency relationship between workstations determined by the path of the automated guided vehicle and the process flow is abstracted as a directed edge. The direction of the edge is from the preceding workstation to the following workstation. The edge weight represents the sensitivity of the completion delay of workstation i to the start time of the following workstation j. The edge weight is calibrated by the regression coefficient of the delay of workstation i and the waiting time of workstation j in the historical production line operation data, or by measuring the offset of the start time of workstation j under different delay amounts through digital twin simulation. The calculation of the cascaded disturbance propagation intensity includes: when any workstation i experiences an anomaly and generates an initial delay time, it propagates along the directed edge to the subsequent workstation j. The delay increment after propagation is obtained by multiplying the upstream delay by the edge weight, subtracting the preset dynamic buffer time of the downstream workstation, and comparing it with zero to take the larger value. A rule was set that if multiple paths lead to the same workstation, the maximum value of the delay increment calculated for each path would be taken as the disturbance to that workstation. According to the above rules, the disturbance propagation proceeds step by step from the source workstation of the anomaly to the downstream of the directed graph until it reaches the workstation with zero out-degree in the directed graph of the production line disturbance propagation. This workstation is the end workstation of the production line. The delay increments on all propagation paths of the end workstation of the production line are accumulated to obtain the cascade disturbance propagation intensity corresponding to workstation i. When no abnormality occurs at workstation i, the cascade disturbance conduction intensity is zero; when multiple workstations experience abnormalities simultaneously, the cascade disturbance conduction intensity is calculated independently for each abnormal workstation as the source.

5. The intelligent management method for a flexible production line of a large fixed-wing UAV according to claim 4, characterized in that: The rescheduling trigger factor is generated by obtaining the cascading disturbance propagation intensity of all workstations. 、 The order urgency correction coefficient and workstation importance weight associated with each workstation are calculated by multiplying the importance weight, urgency correction coefficient and cascading disturbance propagation strength of each workstation, and then summing the products of all workstations. The result is used as the rescheduling trigger factor. The calculation process for the order urgency correction coefficient is as follows: If there are no orders in production at workstation i, the order urgency correction factor is set to zero. If there are orders in production, the difference between the delivery deadline of the order and the current system time is obtained as the remaining available time. The total time required for the assembly process of the order that has not yet been completed at workstation i is obtained as the required processing time. The ratio of the remaining available time to the required processing time is calculated. Obtain the preset maximum ratio threshold in the system and calculate the order urgency correction coefficient. The method for determining the importance weight of each workstation is as follows: calculate the average work-in-process dwell time and process output rate of each workstation in the production line, identify the workstation with the longest average work-in-process dwell time and the lowest process output rate as the bottleneck workstation, assign higher workstation importance weight values ​​to the bottleneck workstation and its immediate predecessor workstation, assign the lowest workstation importance weight values ​​to non-bottleneck workstations with parallel substitution capabilities, and normalize all weights to a sum of 1. The preset trigger condition is: when the rescheduling trigger factor is greater than the adaptive threshold, dynamic rescheduling is initiated; The adaptive threshold calculation process is as follows: Set a time window of length L, which slides forward with the system time with a preset step size. Obtain all historical rescheduling trigger factor samples within the window, calculate the sample mean and sample standard deviation, and add the product of the sensitivity coefficient and the standard deviation to the mean as the threshold. The sensitivity coefficient is preset according to the production line's tolerance to disturbances.

6. The intelligent management method for a flexible production line of a large fixed-wing UAV according to claim 5, characterized in that: In step S5, the assembly datum adaptive compensation parameters are generated in the following way: When rescheduling leads to tooling replacement or product configuration switching, the positioning deviation data of key feature points of the component to be assembled is obtained from the feedback of the online measurement system. The current working condition parameters and positioning deviation data are input into the pre-built compensation prediction model to obtain the initial compensation amount. In the digital twin scenario, the initial compensation amount is loaded, and the verification is performed using a mechanism model that includes assembly stress and geometric interference verification. The root mean square value of the residual of the key point after compensation is calculated. The calculation process of the root mean square value of the residual is as follows: obtain the deviation between the compensated position and the theoretical position of all key feature points, calculate the square of each deviation, obtain the square value, calculate the average value of the square value, and then take the square root. If the root mean square value of the residual at the key point after compensation is less than the preset acceptance threshold, the initial compensation amount is used as the adaptive compensation parameter for the assembly datum. Otherwise, based on the feedback of the verification deviation, the gradient descent algorithm is used to iteratively correct the compensation amount. The root mean square value of the residual at the key point after compensation is used as the loss function, and the compensation amount components of each locator are used as optimization variables. The compensation amount is updated along the negative gradient direction. The update step size is preset according to the locator resolution and is constrained by the physical travel limit of the locator. The iteration continues until the root mean square value of the residual at the key point after compensation is less than the acceptance threshold or the preset maximum number of iterations is reached. Finally, the displacement compensation command value is output.

7. The intelligent management method for a flexible production line of a large fixed-wing UAV according to claim 6, characterized in that: The compensation prediction model is a deviation compensation model based on the configuration adaptive influence coefficient matrix. Its construction and usage process includes: During the production line debugging phase, for each UAV product configuration and tooling combination, virtual assembly simulation is performed through a digital twin scenario. Unit displacement excitation with a preset amplitude is applied to each positioner in sequence. The preset amplitude is selected within the repeatability accuracy range of the positioner. The displacement response of the corresponding key feature points is recorded to form the excitation-response data pair of the combination. For each configuration and tooling combination, virtual simulation is repeated multiple times, and valid data recorded in the actual assembly history is superimposed. The virtual simulation data weight and the actual data weight are weighted and fused. After removing samples whose deviation exceeds the preset abnormal threshold, the unit excitation vectors of each excitation are arranged into an excitation matrix, and the corresponding key point displacement response vectors are arranged into a response matrix R. The weighted least squares method is used to identify the influence coefficient matrix under this combination. The calculation method is to minimize the weighted residual sum of squares. During online prediction, based on the product configuration and tooling used for the current task, the corresponding influence coefficient matrix is ​​retrieved from the influence coefficient matrix library. The key feature point positioning deviation vector obtained from online measurement is used as input, and the compensation vector is used as the variable to be optimized to construct a weighted objective function. The optimal compensation vector is obtained by solving a system that makes the derivative of the weighted objective function with respect to the compensation vector zero. The optimal compensation vector is the initial compensation vector.

8. The intelligent management method for a flexible production line of a large fixed-wing UAV according to claim 7, characterized in that: The optimization objectives of the multi-objective optimization in step S5 include: minimizing the order delay penalty value, minimizing the variance of the working hours load of each workstation, and minimizing the total magnitude of the assembly baseline adaptive compensation parameters; The method for minimizing the order delay penalty value is as follows: for each order, subtract the delivery deadline from its expected completion time. If the result is positive, take the difference; otherwise, take zero. Then multiply by the unit time delay penalty weight corresponding to the order. Finally, sum the above products for all orders. The calculation method for the work hour load variance of each workstation is as follows: obtain the planned total processing time and the rated available time of each workstation after rescheduling, calculate the load rate of each workstation, and then calculate the variance of the load rate of all workstations. The total magnitude of the assembly datum adaptive compensation parameters is the sum of the absolute values ​​of each component in the final compensation parameter vector; The constraints for multi-objective optimization include at least the following: each workstation can only process one product at a time; the product process must follow the preset process sequence; the product start time must not be earlier than the material placement time; and the compensation amount of each positioner does not exceed the upper and lower limits of its physical stroke. The compensation scheme generation process adopts a hierarchical sequence method. Minimizing the order delay penalty P is given as the highest priority. The optimal value of this priority is solved, and all solutions whose target value deviation does not exceed the first preset tolerance are retained to form the current feasible solution set. Within the current feasible solution set, minimizing the variance of the working hours and load of each workstation is searched as the second highest priority. Solutions that satisfy this second-highest value are retained to form the updated feasible solution set. Finally, in the updated feasible solution set where the first two priorities are optimal, minimizing the total amplitude of the assembly benchmark adaptive compensation parameter is given as the final priority. The final scheduling sequence and compensation parameters are then obtained.

Images