L Bettoni
email: l.bettoni@unibs.it

L Mazzoldi
email: laura.mazzoldi@unibs.it

I Ferretti
email: ivan.ferretti@unibs.it

L Zavanella
email: lucio.zavanella@unibs.it

S Zanoni

Integrated Energy Value Analysis: a new approach

Keywords: Energy Flow, Energy Analysis, Energy Value Stream, Energy Audit, Lean and Energy

The improvement of energy efficiency of 20% is one of the three objectives of the EU Directive 20-20-20. To reach this goal, all production processes have to be analysed with reference to their energy consumption so as to identify actions aimed at removing or reducing energy wastes. From the Lean Production framework, the variation of the Value Stream Mapping (VSM), the Energy VSM (EVSM) is selected as a useful tool able to highlight the seven types of wastes identified by Ohno in Toyota Production System and energy wastes. This work aims at proposing a possible modification of the EVSM, which encompasses it within the Energy Audit and the Energy Balance Chart. The goal is to realize a deep energy analysis, highlighting energy wastes, in order to understand which corrective actions must be implemented or which corrections to energy reduction should be considered to reduce energy wastes.

Introduction

The study of energy consumption and costs in energy intensive sectors is of primary importance, because of the increase in energy price and its related environmental concerns. In literature, Energy issue is a wide topic, so it is faced from different point of view. Summarizing some of these, we want to highlight the magnitude of energy subject. Some studies have focused on the assessment models of energy consumptions of the system considered, which have the purpose to collect energy data allowing to highlight highest energy consumption process or equipment and quantitatively evaluate energy losses of an industry process [START_REF] Wei | An empirical analysis of energy efficiency in China's iron and steel sector[END_REF][START_REF] Saidur | End-use energy analysis in the Malaysian industrial sector[END_REF][START_REF] Hong | Energy flow analysis in pulp and paper industry[END_REF]. Other studies deal with the measurement of the on-time energy consumptions, in order to straightaway solve different energy-decision problem, and on the solution of multi-objective energy problem related to energy policies, energy operations and management. As the survey in [START_REF] Zhou | Decision analysis in energy and environmental modeling: An update[END_REF], which states that the Decision Analysis (DA) is largely applied to analyze Energy and Environmental (E&E) issues. Different researches are focalized on energy management system [START_REF] Seyithan | Evaluation of corporate energy management practice of energy intensive industries in Turkey[END_REF] or on technological solutions for reducing energy consumption and improving energy efficiency, in comparison with BATs (Best Available Techniques) standards, even if there is a poor knowledge of BATs [START_REF]Reference document on best available techniques for energy efficiency[END_REF][START_REF] Saygin | Benchmarking the energy use of energy-intensive industries in industrialized and in developing countries[END_REF]. Energy intensive companies needs to identify all type of energy use in their processes, to divide which type of energy use is necessary for the production and which is removable (waste). All energy use, which is not strictly necessary for the production, is a cost that the customer is not willing to pay for. For this reason is necessary to provide to the energy intensive companies a useful tool which allows identifying this type of energy use. Even though there are lot contributions in literature, we observed a lack in the availability of an energy analysis able to bring out the different type of energy use. The objective of the Integrated Energy Analysis is to analyse energy use from a top down approach. This approach allows identifying at first macro-evidences of energy use and waste, and gradually highlighting how energy use and wastes are strictly related to production and logistical issues. This research is focused on tools able to identify and quantify energy consumption, waste, energy value and flow in an energy intensive production process. This study, highlighting energy wastes and anomalies and integrating tools derived from Energy and Lean analysis, proposes a new approach (named Integrated Energy Value Analysis) for the analysis of energy features of a production system. Related to other energy analysis, in this research the type of energy waste is pointed out. This information is important to identify the exact action to improve to reduce or remove it. The reduction of wastes is the main objective of Lean production, which aims at reducing costs, improving process flow and reducing lead times, to meet customer expectations and improve environmental quality and in order to achieve competitiveness. To apply the five tenets of Lean production, it is necessary to remove or at least reduce wastes "Whenever there is a product (or service) for a customer, there is a value stream. The challenge lies in seeing it" [START_REF] Rother | Learning to See: Value Stream Mapping to Add Value and Eliminate Muda[END_REF]. Therefore, to see the value stream of a product through the production process, a tool has been developed: the Value Stream Mapping (VSM). VSM is a suitable tool for waste identification and their subsequent elimination (or almost reduction), thanks to its graphical appeal. Some scientific works, as Erlach [START_REF] Erlach | Energy Efficiency in Manufacturing Using the Energy Value Stream Method for Building an Energy-Efficient Factory[END_REF], face the issue of directly evaluating energy efficiency into Lean approach. Erlach introduces the energy contribution in the VSM including energy indicators such as Energy Intensity and Energy Density into the process box. Moreover, he defines eight design guidelines to derive measures for energy efficiency increase. Bogdanski [START_REF] Bogdanski | An Extended Energy Value Stream Approach Applied on the Electronics Industry[END_REF] proposes an extended energy value stream approach, which introduces three innovative characteristics. The main novel feature is the introduction of the state of the machines, classified as ramp-up, processing and idle. For each state the measurement of power and time are collected, in order to consider dynamic energy profile. The second feature is the consideration of energy consumption of technical building service. The third distinction of the extended EVSM is the consideration of product related issues as technological characteristics and design. Another guideline about energy improvement through Energy Value Stream Analysis is the approach of the "Lean and Energy Toolkit", published by U.S. Environmental Protection Agency [START_REF]The Lean and Energy Toolkit[END_REF]. Efforts on Energy Value Stream Mapping are needed.

The remainder of the paper is organized as follows: in section 2, the Integrated Energy Value Analysis methodology is introduced. In section 2.1, the proposed Energy Audit approach is presented; the Energy Value Stream Mapping is explored in section 2.3 and in section 2.4 the Energy Balance Chart is studied in depth. The Future State is described in section 4. The conclusions and further developments are summarized in section 5.

The Methodology

The Integrated Energy Value Analysis (IEVA) is an integration of tools for the energy consumption assessment of an energy intensive process. The focus is on energy issues, so logistical performances are not encompassed by energy analysis and considerations.

Next to the seven wastes identified by Ohno [START_REF] Womack | Lean Thinking: Banish Waste and Create Wealth in Your Corporation[END_REF], during the years researchers added some other wastes, considering different features of the system investigated and the need and criticism to the system itself, i.e. applying Lean and VSM to the service companies [START_REF] Pavnaskar | Classification scheme for lean manufacturing tools[END_REF]. Power and energy are considered as additional wastes, but in the Integrated Energy Value Analysis, energy is the topic of the research. The IEVA analyses different kinds of energy consumptions (and wastes) related, e.g., to the state of the machines or the synchronization of the resources. The objective of the IEVA is to find the hidden energy wastes existing in the production process in order to remove, if possible, or reduce them by energy scheduling models.

Thanks to the analysis of the energy consumption of the machines, different states corresponding to different power levels are observed. In the production state, the machine produces and the power is consistent with the production effort. In the holding state, the machine doesn't produce, but the item newly worked out has not yet been released (i.e., the power is lower than in the production state, as the machine simply holds the item). In the idle state, the lowest power level, the machine is turned on but it is neither producing any items (i.e. no items are in process) nor holding them: nevertheless, some energy is still consumed due to auxiliary equipment of machines, which must be kept in function.

As for the VSM, the energy consumption can be Value-Added, (VA), Non Value-Added (NVA) and Necessary but NVA (NbNVA).

We identified different definitions of "energy waste". (a) When it is possible to consume less energy to produce the same output (i.e. quantity of produced products). (b) When energy is employed for NVA and NbNVA activities. And finally, (c) when energy is consumed supporting Lean thinking wastes, i.e. energy consumed for defective products, energy consumed for useless transportation, energy consumed for extra processing, defective production.

The IEVA focuses its attention on the first two definitions of energy wastes, i.e. (a) and (b), because the other ones are reduced or removed acting on Lean wastes. Thus, the energy wastes, which the IEVA focuses on, are:  Energy consumption during Idle State (i.e. during the set-up),  Energy consumption during Holding State,  Energy Consumption for Overproduction,  Energy lost for lack of Synchronization between energy-consuming resources and their associated auxiliary services. The objective of the investigation is to create an "as is" analysis in order to verify energy consumptions, energy costs and energy wastes. The goal is to reduce both energy consumptions and costs, optimizing energy management, producing more products and consuming less energy, thus becoming "energy efficient". The IEVA allows the identification of energy wastes in order to find corrective actions able to remove or reduce them. The first step is to implement an Energy Audit, distinguishing processes and services resources and classifying resources as Value Added, Non Value Added and Necessary but Non Value Added. The second steps is to create the Energy Value Stream Mapping with related Energy Lines and study the wastes at each single resource (EVA-ENVA); study the synchronization between resources and related auxiliary services (vertical analysisfocus on the single machine). The last step is to create the Energy Balance Chart in order to compare takt time and energy consumption for those process resources, which are more energy consuming. At the end analyses the EBC in order to define an energy scheduling model suitable for the energy wastes reduction, considering the production time value (under/over takt time). In the following sections, the tools are analysed in detail.

Energy Audit

As first step of the IEVA methodology, the Energy Audit (EA) is performed. This technique uses the Lean organization approach in order to look for wastes in the energy use. The philosophy underlying this approach is based on the understanding of how energy is used in a particular process in order to discriminate what is strictly necessary (Value Added) and what can be improved (reduction or elimination of Non Value Added). The phases of the project can be summarized in (1) Energy Analysis, (2) Management analysis and finally (3) Root cause.

During the first phase the data of energy consumption (i.e., rated power, hours of operation per year, estimate of load and contemporary factors and efficiency) are collected and the consumptions at the resources are divided between electrical and thermal loads, with the aim of reconstructing the energy consumptions and costs, reporting a yearly energy balance. The energy balance allows the user to appreciate the annual energy consumption at each user (resource) and a Pareto analysis can be performed to identify the more expensive utilities in the energy perspective.

The Management analysis allows the management to identify system bottlenecks, by the analysis of production processes and the estimation of lead times, thus establishing the actual production rate of the system. The comparison between this value and the production capacity at each station within the process allows the identification of the non-utilized production capacity at the different stations as compared to the bottleneck. This assessment is very important because energy cost, as many other types of cost, can be divided into a fixed component and a variable component, and such components can be identified thanks to this specific activity. Moreover, the resources of the system are further classified as VA and NVA, by the analysis of the states of the machines (production state, holding state and idle state).

Finally, during the last phase of EA, technological improvements are identified and prioritized thanks to the Root cause analysis, which allows us to verify any possible improvement actions, understanding the motivations that underlie the costs currently incurred by the company.

Energy Value Stream Mapping

The EVSM is structured into three sections, as presented in Fig. 1: the Process of production flow (at the top), the Auxiliary service flows (in the middle) and the Energy lines (in the bottom). The EVSM is configured as a particular modification of the traditional VSM: changes are implemented in order to focus only on energy. In the process box of the EVSM, different values are specified: the energy consumption during production, the energy consumption during holding time and the uptime, i.e. the idle time of the process. All auxiliary services necessary for the production are shown in another process box, under the main process stream, and once again energy consumption and uptime are specified. In the IEVA, instead of visualizing the time line, the Energy Lines is shown, distinguishing an energy value added line, EVA, (for the production state) and energy non value-added lines, ENVA (for the holding and idle states). Using this configuration, two analyses can be carried out, in order to obtain two kinds of outputs: (1) Identification of single resource's wastes: the inefficiencies of the single resource are analysed and (2) Identification of wastes related to the synchronization between process resource and related auxiliary service, which works together. Once the removable inefficiencies have been identified (so ENVA), this tool allows the user to define corrective action plans which enable to create the Future State. Fig. 2. EVSM and Energy Line Fig. 3. shows an example of the energy profile of a resource which changes its state over time.

Fig. 2. EVSM and Energy Line

Fig. 3. Example of a waste for lack of Synchronism

In the EVSM, as well as the traditional VSM, the Future State shows the goals the company aims at in the future, i.e. the energy consumption has to be reduced during the holding and idle states, which are ENVA, with respect to the Current State. The most suitable method that can be applied to reduce holding and idle times has to be investigated using energy scheduling models, the development of which is not in the aim of this work. The second analysis, allowed by the EVSM, includes both kinds of resources: process resource and auxiliary service resource. The aim is to verify whether the auxiliary service resource is synchronized with its associated process resource or not. Often, in companies it is observed that the auxiliary service resource (for example, the fun) is in function even if it is not necessary for the associated process resource (for example the furnace), which is, e.g., in its holding state (or in its idle state). In this case the resources are not synchronized, Fig. , and an energy waste is present.

A possible Future State would be that the auxiliary service is turned off (or set to a minimum power level) when the process resource doesn't need it.

Energy Balance Chart

Another tool, derived from the VSM, is used during the Current State definition step, in order to identify other kinds of wastes and to determine the optimal configuration of the Future State: such a tool is the Balance Chart (BC). In the traditional BC, once defined the Takt Time, production times at all the resources (machines or operators or the production times to produce different products in a single machine) are pointed out, in order to identify the bottleneck of the system. Once located the bottleneck (which identifies a type of waste, as it breaks off the flow of the value at takt time), Lean thinking provides some solutions based on system's features, aiming at balancing the flow of the production mix. Thus, the BC allows the analysis of the Current State situation, highlighting the takt time, the resource production time, which has to be lower than the takt time, and the bottleneck of the system, for which the production time is higher than the takt time. In the energy perspective, there isn't an indicator suitable for quantifying how much energy is needed to define an energy efficient system and, at the same time, to maintain the logistic performances linked to the takt time. Therefore, the objective is to find the optimal (i.e. the minimum) energy consumption, in order to comply with both the takt time and the due date required by the customer. To build up an Energy Balance Chart it is necessary to analyze the energy consumption for the current cycle time, for that process resource that is more energy consuming. In order to find the cycle time that ensures the minimum energy consumption, the Total Energy function has to be minimized. The optimum found is the indicator that quantify how much energy is needed at the resource, maintaining the logistical performance. If it is possible to act on the production times of the machines, the Energy Balance Chart helps to highlight the energy behaviour of machines themselves and to identify the optimal production and holding times to reduce energy wastes. As for the traditional Operation Machine Production Balance Chart, the operation machine production time may be lower or higher than the takt time.

To carry out the analysis, we consider one resource or machine (M1): two different cases related to production times may be defined: (1) production time under takt time and (2) over takt time. The following machine, namely M2, which has to be replenished by M1, on the contrary, it runs at the takt time and there is no intermediate buffer between the two machines. So, in the first case (under takt time) the machine M1, once ended the production of the current item, stays in the holding state till the fol-lowing machine becomes available to receive the item itself, i.e. till the takt time. In the second, M1 finishes the production later than the job completion by machine M2, which has to stay in idle state from the end of its production time (i.e. takt time) till M1 finishes to produce its item. The analysis is carried out in order to find the minimum energy consumption related to machine M1.

We assume that the production rate is variable, and that the energy consumption for production can be described as a convex function of the production time, while energy consumption for holding, on the other hand, decreases linearly when increasing the production time, i.e. it is proportional to the holding time.

Let us introduce the following notation in order to define the mathematical functions of energy consumption for production and holding states: We also assume that the energy consumption for production, Ep, is a function of the production time, as in [START_REF] Tang | A stochastic production planning problem with nonlinear cost[END_REF] by: ,

while the energy consumption for holding, Eh, is defined by: when .

The holding energy consumption is proportional to the holding time and decrease linearly increasing production time.

The Total Energy Consumption is the summation of energy consumption for production, Ep, and energy consumption for holding, Eh, and it is expressed by: .

The energy consumed when the resource runs at takt time is defined as: .

The optimal combination of production time and holding time, which minimizes energy consumption, can be found by setting the first derivative of the total energy consumption function to zero, from which the optimal production time that minimizes energy consumption is derived as:

.

(5) The optimal production time tp * is determined for the resources studied, i.e. the production time that minimizes the Etot, The EBC shows the distribution of energy consumptions for different production time values. Thanks to the EBC, it is possible to compare the current Total Energy Consumption with the minimum Total Energy Consumption, which represents the goal to be reached to become "energy efficient". It has to be noted that the EBC identifies the optimal energy trade-off between energy consumption for production and holding. In this study the logistic costs are not considered, so a trade-off between energy savings and logistic costs should be deeply investigated using additional analyses, i.e. energy scheduling models.

Future State

Once defined the energy consumption of the resources and the associated auxiliary services, through the Current State, it is possible to identify the areas of improvement, both management and technological. The type of corrective actions have to be investigated through different models, as the energy scheduling model. Once the results have been achieved from the Energy Value Stream Mapping and the Energy Balance Chart, thus determining the optimal solution, these tools can be revisited in order to edit and improve the Future State, as well as to highlight the energy goal of the system.

Conclusion

The proposed Integrated Energy Value Analysis (IEVA) is a method that allows the managers to investigate energy wastes at different levels of analysis. The Energy Audit, the Energy Value Stream and the Energy Balance Chart are the tools used to perform the IEVA. For the energy consuming resources, different states are classified (production, holding, idle) and each state is categorized as Energy Value Added or Energy Non Value Added. The evaluation of energy consumption is carried out for each single resource to analyze energy use during different states, and for the synchronisation of both resources and associated auxiliary services to highlight the hidden energy wasted because of the incorrect management of resources. The Energy Balance Chart, which allows to relate takt time and energy consumption, in order to understand the optimal energy scheduling model to apply so as to reduce energy consumption. Possible further developments of the present study are the introduction of energy quantitative and qualitative performance measures, in order to control the energy performances, as well as logistical ones. Finally, the definition of a methodol-ogy aimed at selecting the optimal energy scheduling model to be applied has still to be studied in deep.

Fig. 1 .

 1 Fig. 1. EVSM configuration.

  tt : takt time [time unit]; tp : production time of the resource [time unit]; th : holding time of the resource [time unit], defined as: tp, th = tttp α,β, γ : coefficients of energy consumption for production function; δ : coefficient of energy consumption for holding function; Wp : power required for production [kW]; Wh : power required for holding [kW];

  Fig. . Consequently, the optimal configuration of production time/holding time is defined, too. So the Energy Balance Chart (EBC), Fig. , can be drawn for different energy configurations, in order to compare the Current State of the resource studied and to determine the optimal configuration.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Total Energy Consumption


Swee S Kuik
email: kuik@kaede.cs.kobe-u.ac.jp

Toshiya Kaihara
email: kaihara@kobe-u.ac.jp

Nobutada Fujii
email: nfujii@phoe-nix.kobe-u.ac.jp

Daisuke Kokuryo
email: kokuryo@port.kobe-u.ac.jp

An Integrated Production Planning Model with Obsolescence and Lifecycle Considerations in a Reverse Supply Chain

Keywords: reverse supply chain, remanufacturing, obsolescence

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Manufacturers are now aware that product disposal to landfill may not be a viable solution to cope with the increased virgin material and disposal costs [START_REF] Pahl | Integrating deterioration and lifetime constraints in production and supply chain planning: A survey[END_REF][START_REF] Kuik | Sustainable supply chain for collaborative manufacturing[END_REF][START_REF] Nagalingam | Performance measurement of product returns with recovery for sustainable manufacturing[END_REF]. In supply chain perspectives, the depreciation due to its value loss for remanufactured products has a great influence on the profitability of manufacturers. Rapid improvement and development in science and technology has also resulted with shorter product lifecycle [START_REF] Zhao | Varying Lifecycle Lengths Within a Product Take-Back Portfolio[END_REF]. An appropriate production planning for shorter lifecycle products can help reduce the inventory holding costs and improve an overall supply chain performance. In addition, the remanufactured products could lose their perceived utility within end-of-life (EOL) stage [START_REF] Pahl | Integrating deterioration and lifetime constraints in production and supply chain planning: A survey[END_REF][START_REF] Kuik | Sustainable supply chain for collaborative manufacturing[END_REF][START_REF] Nagalingam | Performance measurement of product returns with recovery for sustainable manufacturing[END_REF]. To address these issues of obsolescence, manufacturers should consider to integrate this critical aspect in the time-period production planning horizon.

In response to the significant changes in environmental rules and regulations, numerous manufacturers are also keen for implementing environmentally conscious manufacturing and product recovery with their supply chain partners for achieving significant amount of cost savings as a whole [START_REF] Kuik | Sustainable supply chain for collaborative manufacturing[END_REF][START_REF] Nagalingam | Performance measurement of product returns with recovery for sustainable manufacturing[END_REF]. Remanufacturing is currently known as one of the compromising alternatives to meet the stringent environmental challenges. By considering the alternative dispositions for handling EOL products, manufacturers can then minimise significant amount of waste disposal to landfill. There are organisations that are currently working on product recovery operations, such as Caterpillar, HP, Xerox, etc [START_REF] Kuik | Sustainable supply chain for collaborative manufacturing[END_REF]. For successful production planning, manufacturers need to consider the planning horizon with the obsolescence constraints that may directly impact on the depreciation value of remanufactured products.

At present, the development of production planning with obsolescence consideration for manufacturers is still at the budding stage [START_REF] Pahl | Integrating deterioration and lifetime constraints in production and supply chain planning: A survey[END_REF][START_REF] Nagalingam | Performance measurement of product returns with recovery for sustainable manufacturing[END_REF]. A major limitation is the complication of various types of remanufactured products and component obsolescence constraints, where it is subject to the depreciation value over a certain time-period planning horizon. If the obsolescence restriction is exceeded, any excessive inventory for obsolete products and components may need to be disposed entirely [START_REF] Pahl | Integrating deterioration and lifetime constraints in production and supply chain planning: A survey[END_REF][START_REF] Kuik | Sustainable supply chain for collaborative manufacturing[END_REF][START_REF] Nagalingam | Performance measurement of product returns with recovery for sustainable manufacturing[END_REF]. However, the total cost associated with the disposal treatment is considered as the financial burdens for manufacturers. The reduction of substantial amount of waste disposal by manufacturers is indeed a need for cost-effective improvement in a reverse supply chain.

In this article, an integrated remanufacturing production planning model is proposed with the considerations of component obsolescence restrictions. This developed model may help manufacturers to minimise overall supply chain costs. Section 2 presents the formulation of mathematical models in remanufacturing production planning problems. To demonstrate the usefulness of the proposed models, Section 3 and 4 present the numerical case examples with results analyses. Finally, the contributions of this study and future work are also discussed.

Mathematical Formulation

To address the obsolescence issues in production planning, the mixed integer programming (MIP) optimisation models are developed to analyse various scenarios of costeffective production planning. In this study, there are two mathematical models with or without the incorporation of obsolescence issues that are presented for comparisons. The first model is developed for an analysis to exclude the component obsolescence restrictions. Meanwhile, the second model is developed for an analysis to include component obsolescence restrictions. These optimisation models aim to minimise an overall associated cost of producing remanufactured products for the next planning horizon of T. The detailed mathematical formulations are discussed in the following sections. There are some key assumptions for the developed models that are considered in this study. The input parameters with continuous, integer and binary variables used for model formulation and evaluation are provided in Table 1 and2 (1) Terms 1-3 are the cost associated with product-component in production and inventory aspects. Term 5 is calculated as the set-up associated costs in t period planning horizon. Terms 6 and 7 represent the cost associated with ordering components. Term 8 is the inventory cost associated with component in t period planning horizon. The model constraints are presented in Eqs. ( 2)-(8).

𝐿 𝑗𝑡 = 𝐿 𝑗,𝑡-1 + 𝑆 𝑗 𝑌 𝑗𝑡 -∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡 𝑖 1 < 𝑡 ≤ 𝑜 𝑗 ∀ 𝑡, 𝑗 (2)
𝐿 𝑗𝑡 = 𝐿 𝑗,𝑡-1 + 𝑆 𝑗 𝑦 𝑗,𝑡-𝑂 𝑗 -∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡 𝑖 𝑡 > 𝑜 𝑗 ∀ 𝑡, 𝑗

𝐼 𝑖𝑡 = 𝐼 𝑖,𝑡-1 + 𝑋 𝑖𝑡 -𝐷 𝑖𝑡 ∀ 𝑡, 𝑗 (3)
Eq. ( 2) is the inventory balance associated with order and bill of materials (BOM). Eqs. ( 3) and ( 4) are the formulation of remanufactured products for the inventory balance to satisfy demand request. The constraints with related capacity, such as production and ordering are shown in Eq. (5). Meanwhile, Eq. ( 6) is formulated as the parameters used in binary form, if the values for 𝑦 𝑗𝑡 , 𝑋 𝑖𝑡 , and 𝑌 𝑗𝑡 , and 𝑝 𝑗𝑡 are more than zero, 𝑃 𝑗𝑡 , 𝑎 𝑖𝑡 , and 𝑝 𝑗𝑡 are also more than zero. These values of 𝑦 𝑗𝑡 , 𝑃 𝑗𝑡 , and 𝑝 𝑗𝑡 can be expressed as one with W as a large positive number. In addition, Eqs. ( 7) and (8) represent the integer and binary decision variables used in the modelling.

Model Two with Component Obsolescence and Lifecycle Constraints

This section presents the formulation of MIP model to evaluate production planning with the restriction of component obsolescence. An objective function is formulated with nine terms as shown in Eq. ( 9), where Eqs. ( 10 There are only two modified terms with this proposed model. Term 8 is modified to include the time-period upon receiving components and the end of obsolete timeperiod. The component disposal treatment cost is calculated as shown in the Term 9. The rest of the Terms 1-7 used are similar with the previous model in Eq. [START_REF] Pahl | Integrating deterioration and lifetime constraints in production and supply chain planning: A survey[END_REF].

𝐼 𝑖𝑡 = 𝐼 𝑖,𝑡-1 + 𝑋 𝑖𝑡 -𝐷 𝑖𝑡 ∀ 𝑡, 𝑗 (10)
∑ 𝐿 𝑗𝑡𝑓 𝑡 𝑓 = 𝐿 𝑗0 + 𝑆 𝑗 𝑌 𝑗𝑡 -∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡 𝑖 𝑡 = 1 ∀ 𝑡, 𝑗, 𝑟 (11)
∑ 𝐿 𝑗𝑡𝑓 𝑡 𝑓 = ∑ 𝐿 𝑗,𝑡-1,𝑓 𝑡-1 𝑓 + 𝑆 𝑗 𝑌 𝑗𝑡 -∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡 𝑖 1 < 𝑡 < 𝑒 𝑗 𝑡 ≥ 𝑒 𝑗 ∀ 𝑡, 𝑗 (12)
Eq. ( 10) is the balance of inventory holding for product-component to satisfy demand request. Eqs ( 11) -( 15) are calculated as each time-period of orders is met with scheduled production planning. Eqs. ( 16) and ( 17) represent the production usage at each time-period with the obsolete components. Eqs. ( 18) -( 20) are the expressions of inventory holding for product-component at each time-period. Eq. ( 21) is derived to estimate excessive components to be discarded at the end of the obsolete time-period. In addition, the capacity constraints in relation to the production and ordering are established as shown in Eq. ( 22). Eq. ( 23) represents the decision variable in binary form. If the values of 𝑦 𝑗𝑡 , 𝑋 𝑖𝑡 , 𝑌 𝑗𝑡 , and 𝑝 𝑗𝑡 are more than zero, these variables of 𝑃 𝑗𝑡 , 𝑎 𝑖𝑡 , and 𝑝 𝑗𝑡 are also more than zero. As a result, these variables of 𝑦 𝑗𝑡 , 𝑃 𝑗𝑡 , and 𝑝 𝑗𝑡 are one, where W is a large positive number.

Finally, Eqs. ( 24) and (25) represent the continuous and binary decision variables used in the modelling. In order to resolve the developed optimization models with obsolete component scenario in Eqs. ( 1) and ( 9), a CPLEX mixed integer optimization solver is used.

Numerical Example

This section discusses a numerical case example for production planning problem with the remanufactured products. In this example, a planning horizon is about 12 periods, which is used as a simulated case example. In this remanufacturing production planning example, there are three types (i.e. PT1, PT2 and PT3) of remanufactured product under examination. Each type of the remanufactured products has four separated components, which is named as CP1, CP2, CP3 and CP4. There are considered as the common components that are used by all types of the remanufactured products. Table 3 and4 show that the parameters used for both remanufactured products and components in the modelling. The obsolescence of time-period estimation for each component is about 6 months for CP1, 3 months for CP2, 9 months for CP3, and 3 months for CP4 respectively.

Using Cplex optimisation solver, the results obtained without obsolescence issue showed that an overall minimum cost for the production planning was approximately $5,263,143 with the linear relaxation of approximately $4,829,131. The calculated value of optimality gap is about 0.0825. However, the results obtained with obsolescence issues revealed that an overall minimum cost for production planning was approximately $6,029,983 with the linear relaxation of approximately $5,528,876. The calculated value of optimality gap is about 0.0906. By comparing both models with or without obsolescence, the obtained results revealed that the total cost for production planning model with obsolescence as derived in Eq. (9) would be much higher than the total cost for production planning model without obsolescence as derived in Eq. ( 1). One of the primary reasons is that the product-component inventory holding and disposal are expressed as the cost associated terms of objective functions. In practice, the manufacturer should consider to maintain inventory stocks with quantities of the obsolete components as low as possible to avoid disposal costs after EOL constraints.

Results and Discussions

The obsolescence issues for remanufacturing production planning has been drawn significant attentions in recent years. To minimise the costs of producing the remanufactured products, there is a need to consider the obsolescence issues for product-component level in a production planning horizon. Figure 1 illustrates a summary of product and component inventory variables with the total associated costs for this case example. This consideration may help resolve the oversimplified production planning problems. Two developed models have been compared with the remanufacturing production planning with and without obsolescence considerations. With three remanufactured product types (N=3) and four separate components (J=4), the proposed MIP model with obsolescence scenario consists of 1,107 variables with binary variables of 132, integer variable of 970 and continuous variable of 5. This model contains 2392 constraints used and non-zero coefficient of 1968.

Concluding Remarks

In this study, the proposed production planning model is developed by integrating the obsolescence issues for remanufactured products within a planning horizon. Comparisons of the obtained results from the developed models with and without obsolescence showed that the total associated costs with obsolescence are generally higher than the total associated costs without obsolescence. In addition, the contribution of this study is twofold. Firstly, the developed model addresses the obsolescence issue for remanufactured products in production planning problems. Secondly, the developed model aims to avoid the oversimplification of cost estimation in a reverse supply chain. For future works, different obsolete time-periods for the remanufactured product-component and variations of the inventory holding costs in production planning will be compared. Furthermore, the models can be extended by considering uncertainty of the consumer's demand requirements.



  respectively. Types of the remanufactured products and components are known;  Demands in the planning horizon are deterministic;  Capacity in production and ordering are deterministic and constant;  Quality aspect is excluded for either products or components;  Considerations of obsolete components are deterministic and known;  Transportation and distribution are beyond the scope of this study.

Fig. 1 .

 1 Fig. 1. A summary of total associated costs for producing 3 product types with 4 components

Table 1 .

 1 Model Parameters

	Symbol	Description
	𝑁	Number of remanufactured product 𝑖, where 𝑖 = 1,2,3, … 𝑁
	𝐽	Number of component 𝑗, where 𝑗 = 1,2,3, … 𝐽
	𝑇	Planning horizon comprising of time period t, where 𝑡 = 1,2,3, … 𝑇
	𝑓	Time period f, upon receiving component j,
	𝑅 𝑖	Unit cost of remanufactured product 𝑖
	𝐶 𝑗	Unit cost of component j
	𝑈 𝑖	Inventory holding unit cost for remanufactured product 𝑖
	𝐻 𝑗	Inventory holding unit cost for component j
	𝐴 𝑖	Fixed set-up cost for remanufactured product 𝑖
	𝐺 𝑗	Fixed ordering cost for component j
	𝑉 𝑖	Production capacity for remanufactured product 𝑖
	𝑀 𝑗	Ordering capacity for component j
	𝛼

𝑗 Disposal treatment unit cost for component j 𝑆 𝑗

Order lot size ordering for component j 𝑂 𝑗

Order lead-time for component j

Table 2 .

 2 Model Variables

	Symbol	Description
	𝑋 𝑖𝑡	Quantity of 𝑖 produced in t period
	𝐼 𝑖𝑡	Quantity of inventory 𝑖 in t period
	𝑦 𝑗𝑡	Quantity of j ordered in t period
	𝐿 𝑗𝑡	Quantity of inventory j in t period
	𝐿 𝑗𝑡𝑓	Quantity of inventory in t period upon receiving j at the end of f period
	𝛽 𝑗	Quantity of disposal treatment j
	𝛾 𝑗𝑟𝑡	Quantity of usage in t period upon receiving j at the end of f period
	𝐷 𝑖𝑡	Demand 𝑖 in t period
	𝐵 𝑗𝑖	Bill of Material (BOM) with j of i
	𝑒 𝑗	Obsolescence with j lifecycle
	𝑎 𝑖𝑡	Fixed setup-up of 𝑖 in t period (binary)
	𝑃 𝑗𝑡	Fixed (scheduled) ordering of 𝑗 in t period (binary)
	𝑝 𝑗𝑡	Fixed ordering of 𝑗 in t period (binary)

Model One without Component Obsolescence and Lifecycle Constraints This

  section presents the formulation of MIP model to evaluate production planning problem by excluding component obsolescence constraint. Using the above parameters and decision variables from Table1 and 2, the objective function of MIP model is formulated into eight terms as shown in Eq. (1).

	𝑀𝑖𝑛 𝑍 = ∑ ∑ 𝑅 𝑖 𝑋 𝑖𝑡 + ∑ ∑ 𝐶 𝑗 𝑦 𝑗𝑡 + ∑ ∑ 𝐶 𝑗 𝑌 𝑗𝑡
	𝑖	𝑡	𝑗	𝑡			𝑗	𝑡
		+ ∑ ∑ 𝑈 𝑖 𝐼 𝑖𝑡	+	∑ ∑ 𝐴 𝑖 𝑎 𝑖𝑡	+	∑ ∑ 𝐺 𝑗 𝑃 𝑗𝑡
		𝑖	𝑡				𝑖	𝑡	𝑗	𝑡
		+ ∑ ∑ 𝐺 𝑗 𝑝 𝑗𝑡 +	∑ ∑ 𝐻 𝑗 𝐿 𝑗𝑡
		𝑗	𝑡				𝑗	𝑡

  𝑋 𝑖𝑡 , 𝑦 𝑗𝑡 , 𝑌 𝑗𝑡 , 𝐼 𝑖𝑡 , 𝐿 𝑗𝑡𝑓 , 𝛾 𝑗𝑓𝑡 , 𝛽 𝑗 ≥ 0 ∀ 𝑡, 𝑖, 𝑗, 𝑓 (24) 𝑎 𝑖𝑡 , 𝑃 𝑗𝑡 , 𝑝 𝑗𝑡 ∈ (0,1)

	𝛽 𝑗 = ∑ 𝛾 𝑗,𝑡,𝑡+1-𝑜 𝑗			∀ 𝑡, 𝑗	(21)
	𝑡=𝑜 𝑗 -1
	𝑋 𝑖𝑡 ≤ 𝑉 𝑖		𝑦 𝑗𝑡 ≤ 𝑀 𝑗			∀ 𝑡, 𝑗	(22)
	𝑦 𝑗𝑡 ≤ 𝑊𝑃 𝑗𝑡		𝑋 𝑖𝑡 ≤ 𝑊𝑎 𝑖𝑡	𝑌 𝑗𝑡 ≤ 𝑊𝑝 𝑗𝑡	∀ 𝑡, 𝑖	(23)
	𝑡		𝑡-1
	∑ 𝐿 𝑗𝑡𝑓	= ∑ 𝐿 𝑗,𝑡-1,𝑓	+ 𝑆 𝑗 𝑌 𝑗𝑡
	𝑓=𝑡+1-𝑒 𝑗	𝑓=𝑡+1-𝑒 𝑗		𝑡 ≥ 𝑒 𝑗	𝑡 ≥ 𝑜 𝑗 ∀ 𝑡, 𝑗	(13)
			-∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡
	𝑡	𝑡-1	𝑖
	∑ 𝐿 𝑗𝑡𝑓	= ∑ 𝐿 𝑗,𝑡-1,𝑓	+ 𝑆 𝑗 𝑦 𝑗𝑡 -∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡	1 < 𝑡 < 𝑒 𝑗	𝑡 ≥ 𝑜 𝑗 ∀ 𝑡, 𝑗	(14)
	𝑓	𝑓			𝑖
	𝑡		𝑡-1
	∑ 𝐿 𝑗𝑡𝑓	= ∑ 𝐿 𝑗,𝑡-1,𝑓	+ 𝑆 𝑗 𝑦 𝑗𝑡
	𝑓=𝑡+1-𝑒 𝑗	𝑓=𝑡+1-𝑒 𝑗		𝑡 ≥ 𝑒 𝑗	𝑡 ≥ 𝑜 𝑗 ∀ 𝑡, 𝑗	(15)
			-∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡
				𝑖
	∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡 = ∑ 𝛾 𝑗𝑓t			1 ≤ 𝑡 < 𝑒 𝑗	∀ 𝑡, 𝑗, 𝑟	(16)
	𝑖		𝑓
	∑ 𝐵 𝑗𝑖 𝑋 𝑖𝑡 = ∑ 𝛾 𝑗𝑓𝑡 𝑓=𝑡+1-𝑒 𝑗 𝑖		𝑒 𝑗 ≤ 𝑡 < 𝑇	∀ 𝑡, 𝑗, 𝑟	(17)
	𝐿 𝑗𝑡𝑓 = 𝑆 𝑗 𝑌 𝑗𝑡 -𝛾 𝑗𝑓𝑡		𝑓 = 𝑡	𝑡 ≤ 𝑜 𝑗	∀ 𝑡, 𝑗, 𝑟	(18)
	𝐿 𝑗𝑡𝑓 = 𝑆 𝑗 𝑦 𝑗𝑡 -𝛾 𝑗𝑓𝑡		𝑓 = 𝑡	𝑡 ≤ 𝑜 𝑗	∀ 𝑡, 𝑗, 𝑟	(19)
	𝐿 𝑗𝑡𝑓 = 𝐿 𝑗,𝑡-1,𝑓 -𝛾 𝑗𝑓𝑡		𝑓 < 𝑡	𝑡 -𝑓 ≤ 𝑜 𝑗	∀ 𝑡, 𝑗, 𝑓	(20)

Table 3 .

 3 Parameters used product i for modelling

		PT1	PT2	PT3
	𝑉 𝑖	300	700	700
	𝑅 𝑖	115	100	230
	𝑈 𝑖	56	48	77
	𝐴 𝑖	530	780	840

Table 4 .

 4 Parameters used component j for modelling

		CP1	CP2	CP3	CP4
	𝑆 𝑗	4500	4500	10000	10000
	𝑀 𝑗	60	30	6	8
	𝐶 𝑗	7500	22000	35000	23000
	𝐿 𝑗𝑡	4	5	5	3
	𝛼 𝑗	6	8	9	6
	𝐺 𝑗	420	540	680	460
	𝑂 𝑗	30000	22000	5000	0


Martin Benter
email: m.benter@tuhh.de

Hermann Lödding
email: loedding@tuhh.de

Analysis of Manual Work with 3D Cameras

Keywords: labor productivity, 3d cameras, primary-secondary analysis 1

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The production costs of manufacturing companies are an important factor when facing international competition. Labor costs make up a vital proportion of the production costs when the production is characterized by a high amount of manual processes and high wages. Labor productivity is the main objective for controlling labor costs. Its analysis allows finding areas of improvement and appropriate measures to increase productivity. Classical work analyses such as the MTM system or the REFA systematics require expert knowledge and induce relatively high costs for the execution of the analysis. Small and medium-sized companies often have neither the knowledge about the methodologies nor the capacity to analyze the labor productivity in depth.

This paper presents a low-cost and low-effort approach to analyze labor productivity with 3D cameras. The developed methodology tracks the movement of employees in assembly operations. The collected data is then analyzed with a methodology that is based on the primary-secondary analysis.

2

Productivity analysis

Analyzing labor productivity

Productivity is defined as the ratio between the output and the input of a system. In manufacturing areas, the output consists of the produced goods and is often measured in units. The input consists of the resources that were necessary to create the output. Resources are for example mechanical and human work. The ratio between the created output and the particular input is the respective partial productivity [START_REF] Sumanth | Productivity engineering and management: Productivity measurement, evaluation, planning, and improvement in manufacturing and service organizations[END_REF].

The labor productivity of a manufacturing area is thus the relation between the output of produced goods and the staff working time as input. As the output of a pro-duction is determined by the market, usually the paid working time is analyzed to investigate labor productivity.

The assessment of actual working times can be performed by dividing the paid working time into work steps and measuring the times for these individual steps [START_REF] Refa | Datenermittlung. Methodenlehre der Betriebsorganisation[END_REF]. Instead of using measured actual times, some methods use predetermined times. Basic idea of these methods is to divide work processes in small work steps to which then times are assigned based on defined influencing factors [START_REF] Refa | Datenermittlung. Methodenlehre der Betriebsorganisation[END_REF]. Known representatives are MTM [3] or the work factor method [START_REF] Quick | Das Work-Factor-System[END_REF].

Primary-secondary analysis

The primary-secondary analysis is a method developed by Lotter [5] to analyze the labor productivity of assembly processes. Basic principle is the distinction between primary processes (PP) and secondary processes (SP). This is similar to the Lean Management classification of value adding and non-value adding times [START_REF] Liker | The Toyota Way: 14 Management Principles from the World's Greatest Manufacturer[END_REF]. According to Lotter, all efforts that add value are primary processes. For assembly operations, the most important examples are joining processes. Efforts that do not add value, such as the transportation of parts, are secondary processes.

Based on that, Lotter defines the efficiency (E) as the ratio between the duration of all the primary processes and the total duration. The efficiency equals one if the work task only consists of primary processes and is smaller than one if secondary processes occur. Fig. 1 shows a graphical representation of the efficiency.

Fig. 1. Graphical representation of the economic efficiency

The abscissa shows the amount of primary processes and the ordinate the secondary processes. Movements that only consist of primary movements are drawn horizontally (movement 2) and movements that do not add value are drawn vertically (movement 3). The graphical addition of the processes results in the total cost vector of the analysis with the slope φ. The reduction of this angle corresponds to an increase in economic efficiency.

If the definition of primary and secondary processes is applied strictly, only joining movements are primary processes, while all other movements do not contribute directly to the customer satisfaction and are therefore secondary processes. However, this definition would imply that the proportion of primary processes would be minimal and the economic efficiency would be of little relevance. Lotter therefore proposes to classify necessary movements up to a defined minimum length as primary processes even if they do not add value.

3

The Microsoft Kinect as a 3D camera

Motion capturing methods recognize, track and digitalize human movements for further processing [START_REF] Kitagawa | MoCap for artists: Workflow and techniques for motion capture[END_REF]. The methods can be distinguished by the used type of sensors. Optical methods usually consist of a transmitter that emits infrared light, which is reflected by the analyzed object to the receiver. Depending on the mode of operation, markers are used as reflectors [START_REF] Kitagawa | MoCap for artists: Workflow and techniques for motion capture[END_REF]. The tracking method of the Microsoft Kinect belongs to the non-marker-based systems [8], which means that the natural reflection of the body is used. This has the advantage that the user does not have to wear any devices, which might disturb him or her from the work.

Microsoft developed the Kinect for the gaming console Xbox as an inexpensive motion control. Microsoft also published a software development kit (SDK), a collection of tools and documentation, for non-commercial applications. The SDK includes the functionality to recognize human persons from the coordinates that the infrared camera collected. This way, the relevant joints of a person and their coordinates can be identified. The Kinect 1 detects 20 body joints while the updated Kinect 2, which was released in 2014 [8], detects 26 points. If joints and corresponding limbs are represented graphically, a skeleton image is created, that supports motion capture.

The Kinect is capable to identify the joints of the body and to record their coordinates. Velocity and acceleration of the joints are not measured but can be calculated from the coordinates. A deficit in applying the Kinect for productivity analyses is that it cannot identify stopping points between two movements. Furthermore, it is necessary to combine these stopping points to stopping areas to identify joining or material areas. Another difficulty in carrying out a primary-secondary analysis is the distinction between primary and secondary processes. Therefore, a productivity analysis, such as the primary-secondary analysis, is still very time consuming with the Microsoft Kinect and only executable by experts.

Data collection

This chapter describes how the recorded data is prepared for productivity analysis. Section 4.1 describes the experimental set-up. Section 4.2 shows which raw data is recorded and what data can be derived directly. Section 4.3 describes how the stopping points can be determined and Section 4.4 describes how stopping areas can be deduced.

Experimental set-up

To demonstrate the modified primary-secondary analysis with the Kinect, we recorded an exemplary assembly task. Fig. 2 shows the experimental set-up. In this task, the recorded worker assembled small plastic tractors for 20 minutes. Each tractor consists of four parts, which are located in four material boxes. The recorded joint was the right hand. The worker used this hand to pick up materials from the boxes one, two and three and the left hand for box four. The distance between Kinect and the worker was three meters.

Fig. 2. Experimental set-up

Coordinate tracking

To carry out the analysis, firstly the joint coordinates are read out. This is done with a rate of 30 frames per second. Thus, the data collection corresponds to a work sampling with a very high number of recordings. The following data is read for each joint:  x-, y-and z-coordinate  the tracking status of the joints: the information if the joint was identified or not.

If the joint is identified correctly, the coordinates can be used for the further analysis. The coordinates of the right hand joint are displayed graphically in Fig. 3 Left. The velocity and the acceleration can be calculated using the following equations.

v x, n = Δx Δt = x n -x n-1 t n -t n-1 (1) a x, n = Δv Δt = v x, n -v x, n-1 t n -t n-1 (2)
𝑥 𝑛 : x-coordinate at recording n 𝑡 𝑛 : time of recording n 𝑣 𝑥,𝑛 : velocity at recording n 𝑎 𝑥,𝑛 : acceleration at recording n Velocities and accelerations for the other directions can be determined in the same way.

Determination of stopping points

The coordinates, velocities and accelerations can be used to determine stopping points. They are defined as times, at which a joint stops between two movements. Since the Kinect does not deliver them automatically, algorithms were designed to determine these times:

1. Velocity falls below a limit: This algorithm checks whether the velocity of a joint falls below a predefined value for a certain period. This option is particularly suitable for processes in which longer stops dominate. This recognition algorithm fails when the performed movements are very fast. 2. Sign change of acceleration: This algorithm checks whether the sign of the acceleration changes from negative to positive over a certain period. This option is suitable for processes, where only short stops are performed, for example button pressing. There are problems with this algorithm if the movement gets slower and then accelerates again, but no stop occurs. Then, a false stop may be detected. This problem can be mitigated by combining this algorithm with the first one (adding a velocity limit). 3. Change of movement direction: This algorithm checks whether the motion vector of a joint changes by a defined angle over a certain period. This option is of great value when most movements are non-stop. This option is less suitable when stops are performed without changing the direction of the movement. The combined use of the three algorithms leads to a good detection of stopping points as shown in Fig. 3 Right, in which the stopping points are marked black.

Cluster analysis to determine stopping areas

During manual activities such as assembly operations, many stopping points occur. The points have to be combined to stopping areas to allow meaningful interpretations. Thus, the next step is a cluster analysis to identify the stopping areas. In this approach, a hierarchical cluster analysis is performed which consists of three steps [START_REF] Backhaus | Multivariate Analysemethoden: Eine anwendungsorientierte Einführung[END_REF]:

y-coordinate x-coordinate y-coordinate x-coordinate
1. Determining the similarities: A variety of distance measures exists to measure the similarity. For scale properties such as the x-, y-and z-coordinates of the stopping points, the Euclidean distance measure is commonly used [START_REF] Backhaus | Multivariate Analysemethoden: Eine anwendungsorientierte Einführung[END_REF]. 2. Selecting a fusion algorithm: In this step, objects with small distances are merged. We used an agglomerative technique that starts with each object representing a single cluster and then merges these clusters gradually. This merging process can be performed with different fusion algorithms. In our case the single linkage method is used [START_REF] Bock | Automatische Klassifikation: Theoret. u. prakt. Methoden z. Gruppierung u. Strukturierung von Daten (Cluster-Analyse[END_REF].

Determining the number of clusters:

There are many different statistical concepts that determine the optimal number of clusters. For this approach, a maximum distance is defined when merging clusters [START_REF] Backhaus | Multivariate Analysemethoden: Eine anwendungsorientierte Einführung[END_REF]. The outcome of the cluster analysis are stopping areas as shown in Fig. 4 Left. One can see in the figure all the breakpoints of the right hand and their allocation to different clusters. Cluster G, F and A represent the boxes one, two and three from the experimental set-up. Primary-secondary analysis using 3D cameras

This chapter shows how the primary-secondary analysis can be performed by using the Microsoft Kinect. The analysis determines, which activities add value and which do not. The approach is divided into three steps. First, the joining point is determined from the stopping areas derived in Section 4.3. Subsequently, non-relevant movements are identified. The last step is the classification of the relevant movements in primary and secondary movements.

In the primary-secondary analysis by Lotter, the user determines the joining point and the primary range is set based on that point [5]. Alternatively, the joining point can be automatically identified. The Cluster analysis results in a number of stopping areas. In most cases, the stopping area with the most stopping points will also be the joining area. To make sure the right area is selected, the user still has the option to choose the joining area manually. In Fig. 4 Left stopping area C is the selected joining point. The middle of this joining area then will be used as the joining point for the further investigations. This option has the advantage that in contrast to the first option, the determination is based on the recorded data. In addition, different joining points can be defined for both hands.

In the next step, an overview of movements that took place between the different stopping areas is created. The user now examines the movements. He decides whether the movements were relevant to the assembly task (relevant movement) or not (nonrelevant movement). The latter movements are henceforward referred to as tertiary movements. The number and duration of these tertiary movements can be easily computed from the captured data and may be an important area for further improvements. In the experimental set-up only movements between the joining area (Cluster C) and the boxes (Cluster A, F, G) are classified as relevant movements.

In the primary-secondary analysis by Lotter the second important step for classification is to determine the nearest material box or stopping area to calculate the minimum level of movement. This minimum, which is necessary for the joining, has to be defined by the user [5]. Two additional options have been developed:  Nearest stopping area: In this option, the nearest stopping area is determined (area F in fig. 4 Left). The distance from the center of this area to the joining point is used as minimum level of movement (movement 2 in Fig. 4 Right). Prerequisite for this option is that the non-relevant movements were identified correctly before.  Determination by using ergonomic aspects: This option is based on ergonomic aspects, which means biological data like arm length is used to determine a movement distance which can be done quickly without harm for the body. An example is the distance that the hand can cover without moving the upper arm and without having harmful or uncomfortable angles for the elbow joint.

For our analysis, we chose the first option. Following the determination of the minimum level of movement, the relevant movements are divided into primary and secondary parts according to the primary-secondary analysis. To visualize the results of that analyzes the average length and duration are calculated for every pair of stopping areas that form a relevant movement. Fig. 4 Right shows the relevant movements and their primary and secondary parts. These movements now can be drawn into a vector diagram as introduced by Lotter (see chapter 2.2). Fig. 5 shows this result.

Fig. 3 .

 3 Fig. 3. Left: Right hand coordinates; Right: Right hand stopping points

Fig. 4 .

 4 Fig. 4. Left: Determined stopping areas; Right: Movements

Fig. 5 .

 5 Fig. 5. Graphical representation of the vector

Summary

This article presents an approach to analyze labor productivity using 3D cameras. The methodology developed is based on the primary-secondary analysis by Lotter. It divides the movements of a worker into primary, secondary and tertiary movements and thus shows potential for productivity improvement. In comparison to the classic primary-secondary analysis the presented approach offers the following advantages:

 Recording of actual data: The presented method does not assess the ideal, but the actual process, thus showing tertiary wastes in addition to the secondary ones.  Semi-Automated: By using the Microsoft Kinect, the effort of recording the processes is significantly reduced.  Higher Accuracy: The determination of the joining point from the stopping areas is based on real data and can be determined separately for both hands.  Visualization: The graphical representation of the stopping areas and movements can help to directly determine incorrectly positioned parts.

Several aspects can further enhance the method. For example, additional body joints can be examined. In addition, the method could be used to reveal ergonomic potential.


K Harmouch

Y Bernard

L Daniel Geeps

Stator/Rotor Interface Analysis for Piezoelectric Motors

Keywords: tribology, sliding speed, pre-compressive force, friction model, interface reaction forces, universal stator

Based on a literature survey, this paper presents the implementation of analytical modeling approaches of the contact between a stator and a rotor in piezoelectric motors. These modeling approaches are validated by comparison to experimental results obtained on a specially designed piezoelectric motor. Two types of measurements are presented: the displacement of contact point, and the friction force. The second measurement is used to evaluate the friction coefficient and its dependence upon several parameters such as sliding velocity and pre-compressive force. Both the friction coefficient and the displacement of contact point are used as input in the modeling. The comparison of the predicted and measured no-load speed for different pre-compression levels shows a close agreement.

Introduction

Due to the growing demand of high performance motors, with minimum mass, minimum size and high service lifetime, scientists and researchers have started to investigate the potentialities of piezoelectric motors. In order to achieve optimal operation of the motor in terms of torque density, efficiency and lifetime, there is a number of key factors to address.

A literature survey shows that contact factors such as material properties, dimension of friction layer and pre-compressive force are crucial for motor performance ( [START_REF] Ch | Ultrasonic motors technologies and applications[END_REF], [START_REF] Nesbitt | Modeling of a Piezoelectric Rotary Ultrasonic Motor[END_REF], [START_REF] Zhu | Contact analysis and mathematical modeling of traveling wave ultrasonic motors[END_REF]).

In a first part, this paper presents the implementation of these key factors into a new model of friction drive for piezoelectric motors. In a second part, a first set of validation experiments is proposed. These experiments allow studying the tribology characteristics of contact materials in order to analyze the modeling input parameters. A comparison between modeling and experimental results is finally performed.

Friction drive model

The friction drive model of a piezoelectric motor allows predicting its operation characteristics once the vibration displacements of contact points and the properties of contact materials are known.

Note that this paper only deals with the development of friction drive models and does not address the identification of contact point displacements.

There are four basic models of friction [START_REF] Liu | Modeling and control of piezoelectric inertia-friction actuators: review and future research directions[END_REF] from which other models can be derived in order to simulate as accurately as possible the contact interface.

The first model (see Fig. 1.a) is the most classical dry friction model. With tangential force F applied, the contact can keep adherence (stick) until F exceeds particular value from which the contact slides (slip). In this last case, F doesn't vary with the relative velocity of contact points which is defined here as the sliding velocity.

The static/dynamic model (see Fig. 1.b) is based on the fact that the sticking forces of contact are higher than the sliding ones. This model is widely used in the modeling of traveling waves ultrasonic motors [START_REF] Ch | Ultrasonic motors technologies and applications[END_REF], [START_REF] Nesbitt | Modeling of a Piezoelectric Rotary Ultrasonic Motor[END_REF]. The viscous friction model (see Fig. 1.c) differs from the previous models by the fact that F evolves linearly with the sliding velocity. This model is often used in the modeling of the lubricated contact for applications such as guiding devices (bearings, sliders). This model has already been used in impact drive motor with lubricated interface [START_REF] Furutani | Effect of lubrication on impact drive mechanism[END_REF] and in superposed modes piezoelectric motors [START_REF] Ch | Ultrasonic motors technologies and applications[END_REF].

The LuGre representation (see Fig. 1.d) is the simplest model that can integrate all the previously mentioned aspects. This model is used in the friction modeling of several motor types, especially, inertial type motors [START_REF] Altpeter | Friction modeling, identification and compensation[END_REF], [START_REF] De | A new model for control of systems with friction[END_REF], [START_REF] Driesen | Concept, modeling and experimental characterization of the modulated friction inertial drive (mfid) locomotion principle: application to mobile microrobots[END_REF].

In addition, some friction interface modeling of piezoelectric motors are based solely on empirical results for the relation between the friction coefficient and operation parameters such as the sliding velocity, vibration amplitude, excitation frequency, normal force, etc. [START_REF] Li | Modeling and analysis of stick-slip motion in a linear piezoelectric ultrasonic motor considering ultrasonic oscillation effect[END_REF], [START_REF] Chowdhury | The effect of amplitude of vibration on the coefficient of friction for different materials[END_REF].

Friction model: implementation and results

The proposed model is represented in Fig. 2: the stator and the rotor (slider) are considered as two rigid blocks between which is located a spring with equivalent stiffness of contact K sr and equivalent damping coefficient C r . Attached to the stator, the spring can slide on the rotor face with friction coefficient µ.

The rotor of mass M r is subject to external load force F r , damping force F c , inertial force due to the rotor acceleration ɣ r , normal load P and friction force f. On the opposite side, the contact points of the stator perform an elliptical motion of normal axis u z and tangential one u x . Due to the stator/rotor friction, the rotor moves at a speed V r .

According to the fundamental principle of dynamics:

𝑓 ̅ + 𝐹 𝑐 ̅ + 𝐹 𝑟 ̅ = 𝑀 𝑟 𝛾 𝑟 ̅ (1) 𝑀 𝑟 𝑉 𝑟 ̇+ 𝐶 𝑟 𝑉 𝑟 ̅ = 𝐹 𝑟 ̅ + µ𝑃𝑠𝑖𝑔𝑛(𝑉 𝑠 ̅ -𝑉 𝑟 ̅ ) (2)
The friction drive model is implemented into Matlab-Simulink (see Fig. 3). Input data are: Fr, u z , u x and the interface material properties (Young modulus and Poisson ratio for the rotor (E r , ν r ) and for the stator (E s , ν s ) and friction coefficient µ). K sr is a function of contact materials elasticity and of the shape and dimensions of the contact surface [START_REF] Johnson | Contact mechanics[END_REF].

As we have seen in the previous section, the choice of the friction model is based on the operation mode of the motor. As a first approximation, µ is considered following the static/dynamic friction model (see Fig. 1.b).

The results are shown in Fig. 4 and5. P 0 is the precompressive force applied to maintain contact between the stator and the rotor. In this paper C r is not taken into consideration. Fig. 4 shows the thrust/ velocity characteristics of the rotor for different friction materials. For the corresponding operation conditions, it can be seen that, by using the Plexiglas, the motor gives the best performance (in terms of thrust and velocity). It must be mentioned that the characteristics of the motor are highly impacted by the operation conditions such as the pre-compression (see Fig. 5), elliptical motion of contact, rotor mass, etc.

Experimentation

The conducted experiment is done on a specially designed piezoelectric motor (see Fig. 6,7)

The universal stator is a triangular shaped structure. Two piezoelectric actuators are placed along two sides of the triangle. It is capable of creating a vibrating elliptical motion at the tip of the stator on which the rotor is frictionally in contact. This kind of displacement is ensured by quadrature phaseshifted excitations of two piezoelectric actuators.

The test bench (see Fig. 6) is designed to achieve measurements such as the tip displacements (u x , u z ) and contact reaction forces (P, f). The interesting variables to test are the friction materials (E r , E s , ν r, ν s, µ), the supply voltage (Amplitude U 0 , Phase φ, pulsation ω) and the pre-compressive force P 0. :

Tip displacement measurement

The tip displacement measurements are performed with a laser vibrometer (see Fig. 6.b) [START_REF]OFV-534 sensor head with OFV-2500 vibrometer controller[END_REF] in both directions 1 and 2, normal to the faces of the tip. Therefore, u x and u z are estimated as being the sum of the projection of u 1 and u 2 , along x and z directions respectively.

Friction coefficient measurement

At the contact interface, a reaction force is created by frictional interaction between the stator and the slider. Let P and f be the projections of this reaction force, respectively, along the normal and tangential directions of the contact surface of the slider. By application of the fundamental principle of dynamics (see Fig. 7):

𝑓 = 𝑑 2𝐻 (𝑃2 -𝑃1) (3) 𝑃 = 𝑃1 + 𝑃2 (4)
Therefore, µ can be defined as follows:

µ = 𝑓 𝑃 𝑠𝑖𝑔𝑛(𝑉 𝑟 -𝑉 𝑠 ) (5)

Measurements results

In order to obtain µ, P 0 is applied using precompressing screws. Afterwards, the slider is forced to move along x axis in both directions at a speed V r . A sample of Plexiglas is tested and the results are shown in Fig. 8. With 𝑉 𝑠 = 0 , the equation ( 5 Attributing the linear variation of μ close to the origin to technical limits, we will consider μ as constant (=0.2) with respect to V r . In fact, due to dimensional gaps in the system, close to the end of the slider stroke (where the sliding velocity is low), the contact is loosened and thus f drops.

Similarly, below certain limits, μ is considered independent of P 0 (see Fig. 9). However, the value of μ drops under excessive amount of P 0 .

In order to verify the no-load speed versus precompression motor characteristics obtained by the proposed modeling approach, the universal stator is run at a frequency of 1 kHz under an alternative supply voltage of 60V (pk-pk). Under these loading conditions, we find: U x =0.9 μm, U z =2.4 μm.

The comparison, in Fig. 10, between modeling and experimental results shows a good approximation of V r0 with respect to P 0 . However, at 15N of P 0 and beyond, the motor reaches its limit of operation and experiences blocking situation. This phenomenon is not included in the modeling approach.

Conclusion

In this paper, we have studied the impact of friction material choice on motor performance. In order to model the frictional interaction between the stator and the rotor of a piezoelectric motor, the contact vibration displacements and the properties of contact materials are needed as input of the model. The friction coefficient variation of a pair of materials in contact has been, experimentally, investigated. We can say that the friction coefficient is unrelated to the sliding speed and the pre-compressing as long as the piezoelectric motor is used in a limited range of velocity and pre-compression levels. The modeling and experimental results for no-load speed with respect to pre-compression have been compared with satisfying agreement.

Fig. 1 :

 1 Fig. 1: Basic friction models: a) Coulomb, b) static/dynamic, c) viscous friction, c)LuGre.

Fig. 2 :

 2 Fig. 2: a) Friction drive representation, b) applied forces summary.

Fig. 3 :

 3 Fig. 3: Friction drive Simulink model.

Fig. 4 :

 4 Fig. 4: F r / V r characteristics for several materials of friction layer.

Fig. 5 :

 5 Fig. 5: F r / V r characteristics for different pre-compressive forces.

  𝑟 ). With an amplitude of 21 mm/s for V r , we can see that μ varies along a quasistraight line for V r close to zero, and, quickly, saturates at a value of 0.2.

Fig. 6 :

 6 Fig. 6: a) Universal stator type motor, b) tip displacement measurement.

Fig. 7 :Fig. 8 :

 78 Fig. 7: Force measurement.

Fig. 8 :

 8 Fig. 8: f/P(V r ) for PMMA.

Fig. 8 :

 8 Fig. 8: f/P(V r ) for PMMA.

Fig. 10 :

 10 Fig. 10: No-load speed comparison between modeling and experimental results for different values of pre-compression for PMMA.

Fig. 9 :

 9 Fig. 9: μ(P0) for PMMA.

Fig. 8 :

 8 Fig. 8: f/P(V r ) for PMMA.


Ioannis Z Emiris
email: emiris@di.uoa.gr

Compact Formulae in Sparse Elimination [Abstract of Invited talk at

Keywords: Sparse elimination, permanent, generating function, resultant matrix, discriminant, matrix representation

It has by now become a standard approach to use the theory of sparse (or toric) elimination, based on the Newton polytope of a polynomial, in order to reveal and exploit the structure of algebraic systems. This talk surveys compact formulae, including older and recent results, in sparse elimination.

We start with root bounds and juxtapose two recent formulae: a generating function of the m-Bézout bound and a closed-form expression for the mixed volume by means of a matrix permanent. For the sparse resultant, a bevy of results have established determinantal or rational formulae for a large class of systems, starting with Macaulay. The discriminant is closely related to the resultant but admits no compact formula except for very simple cases. We offer a new determinantal formula for the discriminant of a sparse multilinear system arising in computing Nash equilibria. We introduce an alternative notion of compact formula, namely the Newton polytope of the unknown polynomial. It is possible to compute it efficiently for sparse resultants, discriminants, as well as the implicit equation of a parameterized variety. This leads us to consider implicit matrix representations of geometric objects.

Introduction

It is today a standard approach to use the theory of sparse (or toric) elimination, based on the Newton polytope of a polynomial, in order to reveal and exploit the structure of algebraic systems. Consider Laurent polynomials in N variables:

f i = k i j=1 c ij x a ij , c ij = 0,
where x = (x 1 , . . . , x N ), x e = i x e i i , and {a i1 , . . . , a ik i } ⊂ Z N is the support of f i . Let Q i be the Newton polytope of f i defined as the convex hull of the support.

A more constrained form of structure corresponds to multihomogeneous systems, which are defined on a product

P n 1 × P n 2 × • • • × P n S (1)
of projective spaces: the variables are partitioned into S subsets, or blocks, so that each equation is homogeneous in each block of n j + 1 homogeneous variables. In contrast, semi-mixed systems are those whose polynomials can be partitioned into subsets, or blocks, with fixed monomial support. In sparse elimination, each subset of polynomials has the same Newton polytope. In the case of semi-mixed multihomogeneous systems, the polynomials in each block are homogeneous of the same degree per block of variables.

An interesting example of a semi-mixed multihomogeneous system appears in game theory and was studied in [START_REF] Emiris | Root counts of semi-mixed systems, and an application to counting Nash equilibria[END_REF][START_REF] Mckelvey | The maximal number of regular totally mixed Nash equilibria[END_REF]. Consider a game of S players, each with m i options or strategies. The j-th player plays a totally mixed strategy when the corresponding strategy is chosen randomly with nonzero probabilities p

(j) 1 , . . . , p (j)
m j . A totally mixed Nash equilibrium (TMNE) is a combination of such strategies so that no player can improve his payoff by unilaterally choosing other (pure or mixed) strategy. We express payoff P j of the j-th player choosing the i-th strategy by the polynomial

k 1 ,...,k j-1 ,k j+1 ,...,k S a (j) ...,k j-1 ,i,k j+1 ,... p (1) k 1 • • • p (j-1) k j-1 p (j+1) k j+1 • • • p (S) k S . (2)
Here a

(j) k 1 ,k 2 ,.
..,k S denotes the predefined payoff of the player, assuming each player chooses the pure option k ℓ ∈ {1, . . . , m ℓ }. The equations imply that the payoff P j does not depend on own strategy i, as long as other players do not change their strategies. Eliminating the P j 's leads to a multilinear system with n j = m j -1 in space [START_REF] Bernstein | The number of roots of a system of equations[END_REF]. Actual probabilities are determined from the normalizing conditions

m j i=1 p (j) i = 1, for j ∈ {1, . . . , S}.

Root bounds

Bounds on the number of common roots of an algebraic system admit different types of compact formulae. We start with classic bounds named after Bézout, tight for (multi)homogeneous systems. The most general such bound is the following m-Bézout bound, stated for dehomogenized system. Theorem 2.1 (m-Bézout bound). Consider a system of N equations in N affine variables, partitioned into S subsets so that the j-th subset includes n j affine variables, and N = n 1 + • • • + n S . Let d ij be the degree of the i-th equation in the j-th variable subset, for i = 1, . . . , N and j = 1, . . . , S. If the system has finitely many complex roots in space (1), the coefficient of

x n 1 1 • • • x n S S in N i=1 (d i1 x 1 + • • • + d iS x S ) (3)
bounds the number of roots. For generic coefficients this bound is tight.

We shall give a generating function for the m-Bézout number by extending results in [START_REF] Mckelvey | The maximal number of regular totally mixed Nash equilibria[END_REF][START_REF] Mclennan | The maximum number of real roots of a multihomogeneous system of polynomial equations[END_REF][START_REF] Vidunas | Counting derangements and Nash equilibria[END_REF]]. MacMahon's Master theorem [START_REF] Macmahon | Combinatory analysis[END_REF] is a powerful combinatorial result.

Theorem 2.2 (MacMahon's Master Theorem). [START_REF] Macmahon | Combinatory analysis[END_REF] Let A = (a ij ) be a complex S × S matrix, x 1 , . . . , x S be formal variables, and V denote the diagonal matrix with nonzero entries x 1 , . . . , x S . The coefficient of

x n 1 1 • • • x n S S in S j=1 (a j1 x 1 + • • • + a jS x S ) n j (4) equals the coefficient of x n 1 1 • • • x n S S in the multivariate Taylor expansion of f (x 1 , . . . , x S ) = 1 det(I -V A) (5)
around (x 1 , . . . , x S ) = (0, 0, . . . , 0).

Let us consider semi-mixed multihomogeneous systems, where the number of equation blocks equals the number S of variable blocks, and we have a varying number of equations per block; more importantly, the degree per block is arbitrary. Generically, this system has a finite number of solutions. MacMahon's Theorem yields a multivariate generating function for the m-Bézout bound of these systems, which could be helpful when one seeks root counts for a family of systems.

Theorem 2.3.

[20] Consider a multihomogeneous system on space (1) of

N = n 1 + • • • + n S equations,
where the equations are partitioned into S subsets of exactly n 1 , . . . , n S equations. We assume that for any i, j ∈ {1, 2, . . . , S}, the polynomial equations in the i-th subset have degree a ij in the variables of the j-th variable subset. Let A be the S × S matrix defined by the a ij 's. Then the m-Bézout bound for the multihomogeneous system equals the coefficient of

x n 1 1 • • • x n S
S in the multivariate Taylor expansion of 1/det(I -V A) around (x 1 , . . . , x S ) = (0, 0, . . . , 0).

For TMNE system (2), we have a ii = 0, and a ij = 1 for i = j. Then

I -V A =      1 -x 1 -x 1 • • • -x 1 -x 2 1 -x 2 • • • -x 2 . . . . . . . . . . . . . . . -x S -x S • • • -x S 1      , and 1/det(I -V A) = 1 1 -σ 2 -2σ 3 -• • • -(S -1)σ S ,
where σ j (for j = 2, . . . , S) denotes the j-th elementary symmetric polynomial in x 1 , . . . , x S . Now we pass to the most general root bound in sparse elimination theory, namely mixed volume: there are efficient algorithms but, in certain cases, the mixed volume computation is reduced to a matrix permanent.

The toric root bound is named after Bernstein, Khovanskii, and Kushnirenko (BKK) [START_REF] Bernstein | The number of roots of a system of equations[END_REF]:

Theorem 2.4. For f 1 , . . . f N ∈ C[x ±1 1 , . . . , x ±1 N ]
with Newton polytopes Q 1 , . . . , Q N , the number of common isolated solutions, multiplicities counted, in the corresponding toric variety, which contains (C * ) N as a dense subset, does not exceed MV(Q 1 , . . ., Q N ), independently of the toric variety's dimension.

It is known that mixed volume is related to the permanent. In fact, by reduction to the latter one shows that mixed volume is #P-complete. For a multihomogeneous system on space (1) of N = n 1 + • • • + n S equations, assume the i-th equation has degree a ij in the j-th variable block. Let A = (a ij ) be an N × N matrix with the columns repeated n j times for each j ∈ {1, 2, . . . , S}. Then the mixed volume equals [START_REF] Mclennan | The maximum number of real roots of a multihomogeneous system of polynomial equations[END_REF]Thm 2] 1

n 1 ! • • • n S ! perm A.
A generalization to arbitrary (non-homogeneous) systems follows. For each block of variables there is a n j -dimensional polytope Γ j , each in a separate complementary space, 1 ≤ j ≤ S. Assume the Newton polytopes are direct products of scalar multiples (by a ij ) of the Γ j 's. We obtain an algebraic system on the product of toric varieties corresponding to Γ j , of dimension n j .

Theorem 2.5. [START_REF] Emiris | Root counts of semi-mixed systems, and an application to counting Nash equilibria[END_REF] Let A be the matrix of a ij 's with columns repeated n j times. If there are N = S j=1 n j equations, then

MV(Q 1 , . . . , Q N ) = perm A S j=1 vol(Γ j ), Q i = S j=1 a ij Γ j .

Resultant formulae

For the sparse resultant, a bevy of results have established matrix and, more particularly, determinantal formulae for a large class of systems. An introductory survey can be found in [START_REF] Morozov | New and old results in resultant theory[END_REF].

Resultants provide efficient ways for studying and solving polynomial systems by means of matrices whose determinant is a non-trivial multiple of the resultant. They are most efficiently expressed by a generically non-singular matrix, whose determinant is the resultant polynomial or, when this is impossible, by also specifying a minor M ′ of resultant matrix M which divides the determinant of the latter so as to yield a rational formula for the resultant R: R = det M / det M ′ .

Macaulay's classical result [START_REF] Macaulay | Some formulae in elimination[END_REF] establishes such rational formula for arbitrary dense systems, and D'Andrea, 100 years later [START_REF]Macaulay-style formulas for the sparse resultant[END_REF], for arbitrary systems in the context of sparse elimination. His recursive construction has been simplified in certain cases [START_REF] Emiris | Single-lifting Macaulay-type formulae of generalized unmixed sparse resultants[END_REF].

For two univariate polynomials there are matrix formulae named after Sylvester and Bézout, whose determinant is equal to the resultant; we refer to them as determinantal formulae. The largest family admitting determinantal formulae is a class of multihomogeneous systems studied in [START_REF] Sturmfels | Multigraded resultants of Sylvester type[END_REF][START_REF] Weyman | Multigraded formulae for multigraded resultants[END_REF], see also [START_REF] Gelfand | Discriminants, Resultants and Multidimensional Determinants[END_REF]Sect.13.2], where, each block of variables contains a single dehomogenized variable or is of linear total degree. These systems are moreover unmixed, in other words all polynomials have the same Newton polytope; they are called multigraded.

The multigraded resultant matrices were made explicit in [START_REF] Dickenstein | Multihomogeneous resultant formulae by means of complexes[END_REF][START_REF] Emiris | Multihomogeneous resultant matrices for systems with scaled support[END_REF], which offered determinantal formulae of Sylvester, Bézout and hybrid types, the latter containing blocks of both Sylvester and Bézout types. In [START_REF] Dickenstein | Multihomogeneous resultant formulae by means of complexes[END_REF], it is shown that there exists a determinantal pure Bézout-type resultant formula if and only if there exists such a Sylvester-type formula. The Bézout-type matrices generalize those identified in [START_REF] Chtcherba | Conditions for exact resultants using the Dixon formulation[END_REF]. In [START_REF] Emiris | Multihomogeneous resultant matrices for systems with scaled support[END_REF] was proven the existence of determinantal formulae for systems whose Newton polytopes are scaled copies of one polytope.

In [START_REF] Dickenstein | Explicit formulas for the multivariate resultant[END_REF], hybrid formulae were proposed in the mixed homogeneous case. The most complete study of hybrid formulae is [START_REF] Eisenbud | Resultants and Chow forms via exterior syzygies[END_REF] where resultants are given as a matrix determinant or a Pfaffian, i.e. a determinant square root, for unmixed systems of 3 dense polynomials, for up to 5 polynomials of degree up to 4, 6, or 8, and for 6 quadratic polynomials.

Some studies have focused on three polynomials. In the unmixed case a determinantal formula is established in [START_REF] Khetan | The resultant of an unmixed bivariate system[END_REF], more direct than [START_REF] Eisenbud | Resultants and Chow forms via exterior syzygies[END_REF] and generalizing the Macaulay-style formula obtained for the dense case in [START_REF] Jouanolou | Formes d'inertie et résultant : Un formulaire[END_REF], where the numerator matrix has one row corresponding to the affine Jacobian. For Newton polygons which are scaled copies of a single one, the smallest Sylvester-type matrices are obtained by using a row containing the coefficients of the toric Jacobian [START_REF] Emiris | Hybrid sparse resultant matrices for bivariate polynomials[END_REF]. In [START_REF] Zhang | Rectangular corner cutting and Sylvester A-resultants[END_REF] they construct determinantal Sylvester-type formulae for unmixed systems whose Newton polygon is a rectangle from which smaller rectangles have been removed at the corners.

Discriminant formulae

Discriminants are crucial in studying well-constrained algebraic systems, the system's zero set, its singularities, and other properties of the system. The discriminant is closely related to the resultant but admits no compact formula except for very simple cases. Below we offer a new determinantal formula for a specific system. Recent work has shed light on the degree of mixed discriminants [START_REF] Cattani | Mixed discriminants[END_REF][START_REF] Dickenstein | Plane mixed discriminants and toric jacobians[END_REF] as a function of the Newton polytopes of the corresponding polynomials. An important open question is to obtain an explicit formula for multilinear systems [START_REF] Sturmfels | The Hurwitz form of a projective variety[END_REF].

In the context of sparse elimination the discriminant is defined for a set of polynomials F i with fixed supports

A i , hence it is denoted ∆ A 1 ,...,An (F 1 , . . . , F n ) or ∆(F 1 , . . . , F n ).
Definition 4.1. The mixed discriminant ∆(F 1 , . . . , F n ) of n polynomials in n variables with fixed supports A 1 , . . . , A n in Z n is the irreducible polynomial (with integer coprime coefficients, defined up to sign) in the coefficients of the F i which vanishes whenever the system F 1 = • • • = F n = 0 has a multiple root (i.e., a root which is not simple) with nonzero coordinates, provided his discriminantal variety is a hypersurface; otherwise, the mixed discriminant is equal to the constant 1.

Clearly, the zero locus of the mixed discriminant is the variety of ill-posed systems, i.e., systems with at least one multiple root.

We focus on systems over P 1 × P 1 × P 1 where each equation is missing one variable, motivated by the study of TMNE in games of 3 players with 2 pure strategies each, see polynomials (2). The bilinear system is as follows:

F 1 : a 0 y 1 z 1 + a 1 y 1 z 0 + a 2 y 0 z 1 + a 4 y 0 x 0 = 0, F 2 : b 0 x 1 z 1 + b 1 x 1 z 0 + b 3 x 0 z 1 + b 4 x 0 z 0 = 0, F 3 : c 0 x 1 y 1 + c 2 x 1 y 0 + c 3 x 0 y 1 + c 4 x 0 y 0 = 0.
The projective coordinates (x 1 : x 0 ), (y 1 : y 0 ), (z 1 : z 0 ) represent the unknown probability ratios for the two options of each player; the coefficients a i equal the differences in the given payoffs (between the two options) of the first player, for the 4 combinations of pure choices of the other two players; similarly, b i , c i equal the differences in the given payoffs (between their two options) of the other two players.

The following theorem gives a 6 × 6 determinantal expression for the discriminant. The formula is of Sylvester-type because the matrix entries are either zero or equal to a polynomial coefficient.

Theorem 4.2. [19] The discriminant ∆(F 1 , F 2 , F 3 ) equals det         0 0 c 0 c 1 b 0 b 1 0 0 c 2 c 4 b 3 b 4 c 0 c 2 0 0 a 0 a 2 c 1 c 4 0 0 a 3 a 4 b 0 b 3 a 0 a 3 0 0 b 1 b 4 a 2 a 4 0 0         .

Sparse implicitization

We argue that in sparse elimination it makes sense to consider an alternative notion of compact formula, namely the Newton polytope of the unknown polynomial. It typically allows computation of the polynomial by interpolation of the coefficients. It is possible to compute it efficiently for sparse resultants [START_REF] Emiris | An oracle-based, output-sensitive algorithm for projections of resultant polytopes[END_REF], discriminants, as well as the implicit equation of a parameterized variety.

This section focuses on the latter case. Given (a superset of) the Newton polytope of the implicit equation, one constructs an interpolation matrix, which allows for various geometric operations to be efficiently computed on this matrix without need of developing the monomial representation of the implicit equation. The representation by interpolation matrices has been developed for parametric plane curves, surfaces, and hypersurfaces in the context of sparse elimination [START_REF] Emiris | Geometric operations using sparse interpolation matrices[END_REF][START_REF] Emiris | Sparse implicitization by interpolation: Characterizing non-exactness, and an application to computing discriminants[END_REF]. It generalizes to space curves and, generally, objects of codimension higher than 1 by means of the Chow form, which is a generalization of resultants to N -variate systems with more than N + 1 polynomials [START_REF] Emiris | Matrix representations by means of interpolation[END_REF].

Another powerful matrix representation of geometric objects is based on syzygies. The corresponding theory, including µ-bases, has been developed for parametric models [START_REF] Busé | Matrix-based implicit representations of rational algebraic curves and applications[END_REF][START_REF] Busé | Implicitization of surfaces in P 3 in the presence of base points[END_REF][START_REF] Sederberg | Implicitization using moving curves and surfaces[END_REF]. The computation of syzygies can exploit the sparse elimination setting [START_REF] Botbol | Implicitization of rational hypersurfaces via linear syzygies: A practical overview[END_REF] and, moreover, can be achieved for point cloud models [START_REF] Emiris | Matrix representations by means of interpolation[END_REF]. Roughly, syzygies of a sufficiently high degree define moving curves or surfaces that "follow" the parametric object.

One forms the matrix expressing these moving curves and/or surfaces, which represents the geometric object because its rank at some point p drops if and only if p belongs to the algebraic closure of the image of the parameterization [START_REF] Busé | Matrix-based implicit representations of rational algebraic curves and applications[END_REF].