Microsoft Word - inner cover.docx HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY (HJIC) formerly (until 2012) Hungarian Journal of Industrial Chemistry The HJIC is an international periodicals that focuses on results of fundamental and applied research in the field of • Biotechnology • Chemical Engineering Science • Chemical Processes • Energetics • Environmental Chemistry • Environmental Engineering & Technology • Industrial Management • Material Science • Mechanical Engineering • Mechatronics • Process & System Engineering • Recycling in the form of original papers, reviews, short communications, and conference proceedings written in English. EDITORIAL BOARD Editor-in Chief: RÓBERT K. SZILÁGYI Department of Chemistry and Biochemistry, Montana State University, Bozeman, MT, U.S.A. Honorary Senior Editor: GÉZA HORVÁTH Department of Chemical Engineering Science, University of Pannonia, Veszprém, Hungary Associate Editors: JÁNOS ABONYI Department of Process Engineering, University of Pannonia, Veszprém, Hungary DEZSŐ BODA Department of Physical Chemistry, University of Pannonia, Veszprém, Hungary NORBERT MISKOLCZI MOL Department of Hydrocarbon and Coal Processing, University of Pannonia, Veszprém, Hungary DÓRA RIPPEL PETHŐ Department of Chemical Engineering Science University of Pannonia, Veszprém, Hungary Editors: LÁSZLÓ BARTHA MOL Department of Hydrocarbon and Coal Processing, University of Pannonia, Veszprém, Hungary KATALIN BÉLAFI-BAKÓ Research Institute of Bioengineering, Membrane Technology and Energetics, University of Pannonia, Veszprém, Hungary PETER CZERMAK Institute of Bioprocess Engineering and Pharmaceutical Technology, Mittelhessen University of Applied Sciences, Giessen, Germany DÉNES FODOR Institute of Mechanical Engineering, University of Pannonia, Veszprém, Hungary MARIA GAVRILESCU Department of Environmental Engineering and Management, Gheorghe Asachi Technical University of Iasi, Romania LÁSZLÓ GUBICZA Research Institute of Bioengineering, Membrane Technology and Energetics, University of Pannonia, Veszprém, Hungary JENŐ HANCSÓK MOL Department of Hydrocarbon and Coal Processing, University of Pannonia, Veszprém, Hungary JIŘÍ KLEMEŠ Centre for Process Integration and Intensification, University of Pannonia, Veszprém, Hungary ZOLTÁN KOVÁCS Department of Management, University of Pannonia, Veszprém, Hungary JÁNOS KRISTÓF Department of Analytical Chemistry, University of Pannonia, Veszprém, Hungary ÁKOS RÉDEY Department of Environmental Engineering and Chemical Technology, University of Pannonia, Veszprém, Hungary ISTVÁN SZALAI Institute of Physics and Mechatronics, University of Pannonia, Veszprém, Hungary FERENC SZEIFERT Department of Process Engineering, University of Pannonia, Veszprém, Hungary JÁNOS SZÉPVÖLGYI Research Centre for Natural Sciences, University of Pannonia, Veszprém, Hungary IMRE TÍMÁR Institute of Mechanical Engineering, University of Pannonia, Veszprém, Hungary GYULA VATAI Department of Food Engineering, Corvinus University of Budapest, Hungary GÁBOR VERESS Federation of Technical and Scientific Societies – MTESZ Budapest, Hungary IBOLYA ZSOLDOS Department of Material Science and Technology, Széchenyi István University, Győr, Hungary EDITORIAL OFFICE: UNIVERSITY OF PANNONIA, P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-298, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert Károly PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 Felelős kiadó: a Pannon Egyetem, Mérnöki Kar dékánja EDITORIAL PREFACE Process Engineering and Factory Automation became nowadays the dominant driving forces behind key industries, such as Oil and Gas Industry, Chemical, Pharmaceutical, Food, Machinery, Electronic Industry and Automobile Manufacturing. The aim of the current issue (Volume 41, Number 1) is to overview the latest developments in industrial process automation and process engineering. The papers presented here were selected from contributions at the Factory Automation Conference held on May 21–22, 2013. The objective of the Conference was to establish an international forum for specialists both from the industry and academia. The meeting provided opportunities for the participants to present and discuss the latest developments and up-to-date applications in automation. This issue represents the entire spectrum of process engineering as follows: • Automation tools and solutions of manufacturing processes • Industrial robots in automation • Automation of quality control • Industrial communication networks • Automation safety • Application fields of CAD/CAM, EDA • Lifecycle management • Production control, production logistics • Monitoring and diagnostics of manufacturing processes • Production design, optimization and simulation We are grateful for the authors for their contributions. We anticipate that the readership of the Hungarian Journal of Industry and Chemistry will be inspired and benefit from the selected good quality papers. We plan to continue the tradition of the former Industrial Automation Conference Series. You are welcome to participate at the International Automation Congress (http://www.iac2014.hu) in Budapest, October 29–31, 2014! JÁNOS ABONYI AND ROZÁLIA PIGLER-LAKNER University of Pannonia, Veszprém, HUNGARY Guest Editors Hung. J. Ind. & Chem. 2013, Vol. 41(1), pp. 1–82 Table of Contents Vehicle Routing Approach for Lean Shop-Floor Logistics DAVID GYULAI, LÁSZLÓ MONOSTORI ............................................................................................... 1–6 Design of a New Chemical Injection Pump System for Gas Hydrate Inhibition CSABA VÖRÖS, VIKTOR FÜVESI, ÁKOS PINTÉR ............................................................................... 7–10 Ionic Liquids in Advanced Energy Storage Cells ATTILA GÖLLEI, ATTILA MAGYAR ................................................................................................ 11–15 Monitoring and Diagnosis of Manufacturing Systems using Timed Coloured Petri Nets ADRIEN LEITOLD, BRIGITTA MÁRCZI, ANNA IBOLYA PÓZNA, MIKLÓS GERZSON ......................... 17–26 Immediate Event-Aware Model and Algorithm of a General Scheduler TIBOR DULAI, ÁGNES WERNER-STARK, KATALIN M. HANGOS .................................................... 27–34 Estimation of Parameters of an Extremely Low Fuel Consumption Internal Combustion Engine-Based Megameter-III Vehicle ISTVÁN PINTÉR, MIHÁLY BAGÁNY ................................................................................................ 35–40 Influence of Can Flatness on Heat Dissipation of Aluminium Electrolytic Capacitor LÁSZLÓ KOVÁCS, LÁSZLÓ GÁL, DÉNES FODOR ............................................................................ 41–46 Fine Tuning of Automated Assembly Machines using Video Analysis ALÍZ KATONA, VIKTOR KOVÁCS, DÓRA TASNER, ZOLTÁN KOVÁCS ............................................ 47–50 Quasi-Polynomial Control of a Synchronous Generator ATTILA MAGYAR, ATTILA FODOR ................................................................................................ 51–57 A Structured Model Based Diagnosis Method for Discrete Dynamic Processes using Event Sequences ATTILA TÓTH, KATALIN M. HANGOS, ÁGNES WERNER-STARK .................................................... 59–64 Role of Steady State Data Reconciliation in Process Model Development BARBARA FARSANG, SÁNDOR NÉMETH, JÁNOS ABONYI .............................................................. 65–75 Statistical Quality Control Based Performance Evaluation of Online Analyzers TIBOR KULCSÁR, JÁNOS ABONYI .................................................................................................. 77–82 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 1–6 (2013) VEHICLE ROUTING APPROACH FOR LEAN SHOP-FLOOR LOGISTICS DAVID GYULAI1,2! AND LÁSZLÓ MONOSTORI1,2 1 Fraunhofer Project Centre for Production Management and Informatics, Computer and Automation Research Institute (SZTAKI), Kende str. 13–17, Budapest, H-1111, HUNGARY 2 Department of Manufacturing Science and Technology, Budapest University of Technology and Economics, Egry J. str. 1, Budapest, H-1111 HUNGARY ! E-mail: david.gyulai@sztaki.mta.hu In order to satisfy the material supply needs of large scale shop-floors and production systems, various logistics solutions are applied. In lean manufacturing enterprises, the material supply is pulled by the demands of manufacturing/assembly processes; therefore, a milkrun service is often applied to support the production without glitches. The milkrun logistics planning is a special case for vehicle routing problem (VRP), and requires effective approach to solution in order to satisfy various constraints, and minimize the cost of service. This study gives an overview about lean logistics as well as the most efficient VRP solver algorithms. Furthermore, a novel initial solution with generation heuristics is proposed, which is specially focused on flexible milkrun planning. In order to demonstrate the capabilities of the solution, a software environment is developed as a demonstration that focuses on the main industrial requirements of logistics planning like effective layout definition, quick response of the delivery service and effective order handling. Keywords: vehicle routing problem, milkrun, local search, factory logistics Introduction The material provision of large- and medium-scale production systems requires prudent planning and control, since it influences the performance of the production and affects the order management directly and storage assignment [1]. In order to manage effectively the material provision of the production, practical transport logistics and distribution tasks are usually formulated as vehicle routing problems, whose objective is to obtain a minimum-cost route plan serving a set of customers with known demands [2]. In state-of-the art manufacturing systems, lean processes are successfully applied to eliminate the wastes on all the levels of production, and they raise the overall efficiency of the system by reducing the non- value adding activities in the process chain. The shop- floor logistics processes have to be adapted to the production system; therefore the waste reduction in the material provision also has a key role in balancing the workload and ensuring a smooth running of production. A well-operating lean logistics system keeps the inventories and operational costs on the most cost efficient levels, via by providing the materials what that are needed, exactly when they are needed. For satisfying these requirements, concepts of the lean logistics can be applied that are primarily aimed at controlling the transport services duly to the pull strategy. The paper is focused on the ‘milkrun’ service, which is defined as a manually operated, cyclic transport system delivering raw materials (or sub-assemblies) and disposing empties based on consumption using a fixed route and time schedule [3]. According to this characterization, milkrun scheduling is considered as a special vehicle routing problem with time windows, a limited capacity, and number of vehicles. In this paper a novel, multi-level planning approach is proposed, which combines the advantages of existing algorithms as for example local search for vehicle routing. First, an overview of the state-of-the-art factory logistics is given, considering the recent shop-floor characteristics and planning demands. Then the main solution approaches and algorithms are introduced, as well as the critical points of milkrun planning. In section 4, the proposed layout representation hierarchy and the novel initial solution generation heuristics method are described. The latter focuses on the real industrial criteria and give a feasible milkrun schedule. In the last section, an industrial case-study and some numerical results present the efficiency of the milkrun planning method developed. Shop-floor logistics The shop floor logistics are part of the production logistics in the supply chains of companies. While logistics in general is about supplying a customer with 2 necessary goods in time, with the goal of efficiency, and the lowest possible costs, production logistics is the task of bringing necessary material to workstations and machines and transferring it to the next station. In the realm of logistics planning methodologies the lean approach, originating in the Toyota Production System, of overall waste reduction has become a major planning paradigm, constituting the lean logistics approach [4]. Principles of lean logistics Of the seven waste categories, eliminating unnecessary transports, waiting times, superfluous movement and excess inventory are the most relevant in lean logistics. KLUG identifies the following characteristics of lean logistics [5]: synchronicity and a clocked material flow, flow orientation, standardization, pull principle, stability, integration, as well as perfection. Regarding the pull principle, as the most relevant feature, material is only transferred to a workstation if it is about to process the material thus reducing inventory levels by reducing shop floor stock to safety buffers only. Beside the above principles of lean logistics, a relatively novel issue in material provision of production is ‘green logistics’ that is generally about establishing an ecologically sustainable transportation system considering factors like energy consumption and air- pollution [7]. While in the general external logistics context theses planning principles are reflected in transport concepts such as the just-in-time delivery, there are also corresponding design concepts in production and shop floor logics. For the supplier/source side, establishing zones for decentralized incoming goods and supermarkets in the vicinity of assembly stations is a means of reducing transport lengths. On the receiving side, standardized shelves and reusable containers, replenishment schemes like ‘kanban’ are applied. For the transport system between supply and workstations, ‘tugger trains’ (Fig.1) with fixed time intervals, standardized capacity, and transport lot sizes are preferred over forklift transport in the context of lean logistics; the goal of, which is the reduction of inventory and avoiding erratic peaks in production known as the bullwhip effect [8]. This transport concept of scheduled tugger trains in the context of a lean-oriented logistics system is usually combined with the milkrun transport concept that is in the focus of this paper. Especially in production environments with a high product complexity, variety amid a complex and changing material flow, planning a system of milkrun transports is a complex planning task. Thus, an efficient system of milkruns could significantly in-crease the efficiency of the overall production logistics. Milkrun as a vehicle routing problem Vehicle routing (VRP) is a challenging combinatorial optimization problem, introduced first by DANTZIG and RAMSER in 1959 [9]. The general VRP is known as the Capacitated Vehicle Routing Problem (CVRP), were the aim is to satisfy the needs of all the customers at different locations by having a given number of vehicles with capacity constraints [11]. When milkrun is applied as a material providing system, the time constraint of the runs has also to be added to the general problem (VRP with Time Windows – VRPTW). The general objective is to minimize the costs either by minimizing the total distance travelled, or minimizing the number of vehicles applied. A specific type of VRP is the Pickup and Delivery Problem (PDP) where the vehicles not only provide the locations with materials, but also pick up materials at the stations and deliver them to another one. The milkrun problem is also a PDP, with time windows and cyclic service [12]. Several approaches from the operations science are used to solve VRPs, as for example constraint programming that manages flexibly the various specific constraints of the VRP [14]. Local search is one of the fundamental approaches to find solutions for hard computational – including constraint satisfaction – problems. The basic idea behind local search is to start with a randomly or heuristically generated candidate solution, and to iteratively improve that by means of defined functions [15]. Generating the initial solution is always a crucial point when applying local search, since it significantly affects the running time of the algorithm. Sophisticated heuristics for the VRP are the two-phased ones, which decompose the problem into a clustering- and a routing problem, with possible feedback loops between the two stages. Cluster-First Route-Second (CFRS) algorithms perform a single clustering of the vertex set and then they determine a vehicle route on each cluster. The best- known CFRS algorithm is the FISHER-JAIKUMAR algorithm, which solves the General Assignment Problem (GAP) to form the clusters [17]. The local search starts from the initial solution and subsequently moves from the present solution to a neighbouring solution in the search space where each solution has only a relatively small number of feasible neighbour solutions and each movement is determined by neighbouring operators [20]. When applying local search, the combination of the neighbourhood operators produces the next local optimum solution (Fig.1). Although the most general forms of the VRP can be solved effectively with the above approaches, the milkrun planning problem is hard to interpret as a simple graph-search problem, and the general initial solution heuristics are difficult to transform for this problem. The mapping process of the factory layout also requires special processing; therefore general solution methods cannot be applied easily for such a special problem. 3 Problem formulation Hereby, we provide a definition to the milkrun planning problem. Consider a layout of the shop-floor, which is defined by a set of routes and stations. The routes can be either one- or two-way ones, and each route endpoint is defined as a “routenode”. Along a specific route, there can be more stations (Fig.2). The problem also considers demands that belong to the stations and defined by the amount of the transported goods in standard units and the required cycle-times of the transportation. The goods are transported by vehicles (practically by tugger trains), which are defined by their capacity and average speed. In order to plan feasible milkruns, real-world constraints are considered, such as the capacity constraints of the vehicles, which limit the maximum number of the transported goods. Another constraint is the time limit of the plan, which means the total time consumption of a milkrun plan cannot exceed the cycle- time of the demands. The loading and unloading points can be reached from both sides of the tugger train, and the train can approach the station from both endpoints of the route that the station belongs to. A milkrun plan is built-up by paths, where each path is given as the list of the visited stations and the list of the routes that the vehicle passes along. The plan does not partition the demands, and a visit of a particular station only occurs in a specific path. A milkrun plan is characterized by its total time consumption, which is required by the vehicle(s) to perform the plan. A milkrun plan is considered to be better than the other one only if its total time consumption is smaller, while it satisfies all the demands. The goal of milkrun planning is to minimize the number of the required vehicles via minimizing the time consumption of the plan. The proposed planning process Layout representation for vehicle routing As it was mentioned in the previous section, solving the general VRP is not enough to plan feasible milkrun schedule, since most of the routing algorithms does not calculate with the physical constraints determined by the shop-floor. Most of the general solution approaches consider the problem as a graph-search, where the stations are represented by vertices and the routes are represented by edges [21]. In many cases, this representation does not support effectively the planning processes with feasible results, unless the application of the following constraints. Frequently, when representing a shop-floor by directed graphs, it turns out that the structure of the layout results some direct and indirect routing constraints. A direct case means that the additional equations can be formalized immediately when constructing the graph, such as the asymmetric edge formulas, while the others require further consideration. Indirect constraints are usually implied by the narrow and one-way corridors, which limit the abilities of the vehicles from going to a station directly from another, even if a path connects them. These constraints are implied by the moving abilities of the common applied tugger train, since it cannot turn in the middle of an aisle, and in many of the cases it have to take a detour before visiting the next station (Fig.2). In order to handle the shop-floor constraints effectively and be able to avoid impossible movements of the vehicle, a novel hierarchical layout representation is proposed. The layout is defined by three main different classes, each having their own specific attributes. This structure can be handled dynamically during the path calculations, and the asymmetric nature of the distance matrix is represented together with the limitations of vehicles movements. Fig.3 shows the hierarchy of the layout, defined by “routenodes”, “routes” and “stations”. To calculate the distance matrix of the layout, storing the shortest paths between the stations, DIJKSTRA’S algorithm was applied that solves the single-source shortest path problem in logarithmic running time [22]. In this case, the input of the algorithm is a directed graph where the vertices are the set of the “routenodes” and the edges are the routes. This graph representation is able to handle all the nodes of the layout with their connections, and by applying DIJKSTRA’S algorithm the asymmetric distance matrix can be calculated. Figure 2: Layout elements and vehicle movement Figure 1: Inter- and intra-route relocate operators 4 Initial solution generation In order to plan feasible milkrun cycles with the lowest operational costs, the shop-floor constraints have to be taken into consideration. To handle the constraints effectively, the routing algorithm needs to select the stations to be visited dynamically. Thus, a so-called triplet solution is applied, which helps to avoid impossible movements of the vehicle. A triplet is consisted of three identical nodes: three current, the previous and the next position of the vehicle. The next position is calculated dynamically when determining a particular path of the vehicle. The milkrun planning method strongly focuses on the industrial requirements of vehicle routing; therefore, a novel initial solution generation heuristics was implemented. The goal is to generate a solution that is as close to the criterion as possible. The proposed heuristics can be classified as a CFRS type rule, and differs from the FISHER-JAIKUMAR algorithm by the cluster generation method, since it does not require solving time-consuming GAP, but defines the clusters based on practical reasons [18]. To generate a feasible initial solution, the algorithm calculates a path to each station applying DIJKSTRA’S shortest path algorithm. Then a greedy search algorithm is applied to calculate the next node, which must be visited by the vehicle before returning to the depot, so as to avoid violating the turnaround constraints of the vehicle. The greedy search iterates forward the nodes applying a best-first search strategy [23], and finds the first node from, which the shortest path to the depot does not contain the previous node. The “roundtrip” planning method is applied for each station, and detects all the feasible tours on the shop floor. In order to determine the set of paths for the initial solution, the paths with most visits are selected one by one while the set of unsatisfied demands is not void. Performing the services required by these cycles provides a feasible and acceptable initial solution for the planning algorithm. The pseudo-code of the proposed algorithm is the following: Vv∈ set of stations { }Vv|v:p ∈= path }p{P set of paths [ ] V)s,r(|dD rs ∈= distance matrix }p{I initial solution iq demands of the ith station }q{Q i= set of demands ct cycle time c capacity of the vehicle initialize V, I, Q, C{v}, c, tc foreach i in V calculate D0i foreach v in D0i add v to Ci calculate Di0 foreach v in Di0 add v to Ci align C [length(Ci+1) < length (Ci)] while |Q| > 0 while time(C1) > tc or load(C1) > c remove rand v from C1 (v≠i) foreach v in C1 delete qv from Q add C1 to I delete C1 from C return I Local search strategy and neighbourhood functions In order to improve the initially generated solution, neighbourhood functions are applied. To be able to reach all the points of the search space, both intra- and inter-route operators are necessary to calculate with, although in this case the effect of the inter-route functions is more significant since the initial solution heuristics could provide paths with stations in a rather good sequence. The algorithm performs the local search to minimize the paths, and the cycle-time of them. Then, a simple bin-packing problem is considered, to allocate the paths to a minimal number of vehicles. To solve the problem, first fit heuristics was applied, which is a proven 11/9 OPT solution algorithm [24]. Case study and test results The milkrun planner software, which includes both the user interface and the solver algorithms were implemented in a .NET environment using C++ language. Within the performance evaluation of the implemented software, a real production environment was modelled in order to analyse its capabilities. The purpose of the implemented application is to offer user- friendly and efficient milkrun planning environment; therefore a graphical, point-and-click layout definition interface and XML communication have been applied (Fig.4). Figure 3: Hierarchy of the layout representation 5 First the input parameters are processed and stored in the memory, and DIJKSTRA’S algorithm is applied to calculate the shortest path between each station. The distance matrix is built up from vectors, each having two distance elements. This data structure is required by the applied triplet-based path planning method, which selects dynamically the proper distance parameter from the vector, based on the location of the previously visited station. The initial solution generation uses the distance matrix to create the tours, based on the station loops detected on the layout. The initial solution contains only feasible paths that satisfy the time constraints as well as the capacity constraints of the actual vehicle. To avoid overloaded vehicles, a pre-check method is performed: the vehicle is loaded by all the raw materials in the depot, and then all the other loading and unloading processes are calculated station-by-station to simulate the real milkrun process. In order to calculate with time constraints, both the travel-time (with constant average speed) and the loading/unloading time (time/item) are considered. The purpose of the generated initial paths is to satisfy all the demands by the lowest possible number of feasible cycles. The local search algorithm takes the initial paths and reduces their total time consumption by applying neighbourhood functions where it is possible. The algorithm iterates through all the combinations of the stations, and detects the possible insertions from one path to another. As local search iterates through the paths in a random sequence, different test runs can result different milkrun plans. In order to approximate the global optimum solution, the implemented application performs the local search several times and selects the most appropriate schedule from the generated solutions. The capabilities of the algorithm were tested in a real production environment. This was a large-scale automotive production system with 2 factory halls, 67 stations and 96 routenodes. The milkrun plan has to satisfy 193 various demands in 60 minutes cycle-time. The analysed system requires efficient planning process, since the total length of the routes is over 2.5 km that affects critically the total time of the milkrun schedule. The application of inefficient milkrun schedule results low-utilized vehicles and high number of cycles. In this case, the goal of the test was to minimize the number of vehicles required to perform the services via minimizing the time of the cycles. The proposed initial solution heuristically generated eight feasible paths with 20,400 seconds total time. First, the algorithm tries to optimize the plan by applying one vehicle, and increases the number of the vehicles only if the improvement steps cannot decrease the time requirement of the total plan under the cycle-time (3600 seconds). By this way, the test run could optimize the plan for two vehicles in 177 iterations, and generated a schedule, which time requirement is only 57 minutes (Fig.5). Conclusions After an overview about the state-of-the art shop-floor logistics, and the general formulations of the vehicle routing problem, a solution was proposed that uses a novel layout representation and initial solution generation heuristics to solve the milkrun planning problem. In order to demonstrate the capabilities of the Figure 4: The user interface of the milkrun planner application Figure 5: Test results of the case study 6 solution, a software prototype was developed and tested on real-life industrial data. Future work includes the extension of the model with demand partitioning so as to increase capacity utilization of the transportation vehicles. Moreover, handling inhomogeneous demand types (physical aspects of material handling) requires further constraints to be included in the model. Secondly, future work will be dedicated for making the solution available in practical industrial applications. In an on going research project the milkrun planning algorithm will be implemented in a factory- and a logistics planning application currently under implementation. Acknowledgements The research has been partially supported by the National Office for Research and Technology (NKTH) grant "Digital, real-time enterprises and networks" (OMFB-01638/2009). REFERENCES [1] KOVÁCS A.: Optimizing the storage assignment in a warehouse served by milkrun logistics, International Journal of Production Economics, 2011, 133, 312–318 [2] LAU H., SIM M., TEO K.: Vehicle routing problem with time windows and a limited number of vehicles, European Journal of Operational Research, 2003, 148(3) 559–569 [3] DROSTE M., DEUSE J.: A planning approach for in- plant milk run processes to optimize material provision in assembly systems, Proc. 4th CIRP International Conference on Changeable, Agile, Reconfigurable and Virtual Production, 2011, 605–610 [4] ARNOLD D.: Handbook of logistics, 3rd Edition, Springer, Heidelberg, Germany, 2008 (in German) [5] WOMACK J.P., JONES D.T.: Lean Thinking, Free Press, 2003 [6] KLUG F.: Logistics management in automotive industry, Springer, 2010 (in German) [7] LIN C.: Survey of green vehicle routing problem: past and future trends, Expert Systems with Applications, 2014, 41, 1118–1138. [8] BICHLER K., KROHN R.: Gabler kompaktlexikon logistik, 2nd Edition, Springer, Wiesbaden, Germany, 2011 (in German) [9] SAVELSBERGH M.W.P., SOL M.: The general pickup and delivery problem, Transportation Science, 1995, 29(1), 17–29 [10] DANTZIG G.B., RAMSER J.H.: The truck dispatching problem, Management Science, 1959, 6(1), 80–91 [11] PILLAC V.: A review of dynamic vehicle routing problems, European Journal of Operational Research, 2013, 255(1), 1–11 [12] LU Q., DESSOUKY M.: An exact algorithm for the multiple vehicle pickup and delivery problem, Transportation Science, 2012, 38, 503–514 [13] BERBEGLIA G.: Feasibility of the pickup and delivery problem with fixed partial routes: a complexity analysis, Transportation Science, 2012, 46, 359–373 [14] RUSSEL R.: A constraint programming approach to designing a newspaper distribution system, International Journal of Production Economics, 2013, 145, 32–138 [15] ROSSI F., BEEK VAN P., WALSH T.: Handbook of Constraint Programming, Elsevier, 2006 [16] CACCHIANI V., HEMMELAMAYR V.C., TRICOIRE F.: A set-covering based heuristic algorithm for the periodic vehicle routing problem, Discrete Applied Mathematics, 2014, 163, 53–64 [17] VIDAL T., CRAINIC T.G., GENDREAU M., PRINS C.: Heuristics for multi-attribute vehicle routing problems: A survey and synthesis, European Journal of Operational Research, 2013, 231, 1–21 [18] The VRP Web: http://neo.lcc.uma.es/vrp, (last accessed: May 25, 2014) [19] HASHIMOTO H.: Studies on local search-based approaches for vehicle routing problems, PhD Thesis, Kyoto University, 2008 [20] CARIĆ T.: A modelling and optimization framework for real-world vehicle routing problems, Vehicle Routing Problem, Hrvoje Gold, 2008 [21] MATAI R., SINGH S., MITTAL M.L.: Traveling salesman problem: an overview of applications, formulations, and solution approaches, Traveling Salesman Problem, Theory and Applications, InTech, 2010 [22] WILT C.M., THAYER J.T., RUML W.: A Comparison of greedy search algorithms, Proc. 3rd Annual Symposium on Combinatorial Search, 2010, 129–136 [23] RUSSELL S.J., NORVIG P.: Artificial intelligence: a modern approach, 2nd Edition, Prentice Hall, 2002 [24] YUE M.: A simple proof of the inequality FFD (L) ≤ 11/9 OPT (L) + 1, ∀L for the FFD bin-packing algorithm, Acta Mathematicae Applicatae Sinica, 1991, 7(4), 321–331 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 7–10 (2013) DESIGN OF A NEW CHEMICAL INJECTION PUMP SYSTEM FOR GAS HYDRATE INHIBITION CSABA VÖRÖS1!, VIKTOR FÜVESI1, AND ÁKOS PINTÉR2 1,2Research Institute of Applied Earth Sciences, University of Miskolc, Miskolc-Egyetemváros, POB 2, H-3515 HUNGARY 1Department of Research Instrumentation and Informatics, Miskolc-Egyetemváros POB 2, H-3515, HUNGARY 2Department of Reservoir and Mining Chemical Engineering, Miskolc-Egyetemváros POB 2, H-3515, HUNGARY !E-mail: voros@afki.hu In this study the main feature of chemical injection systems used in gas industry are detailed. The two main energy sources of these systems are air and electric power source. The general structure and main properties of these injection systems are discussed. Examples for a pneumatic commercial and a newly developed, electric system are compared. Keywords: chemical injection system, hydrate inhibition, solar energy, pneumatic system Introduction During the production of gas, a major problem is the formation of hydrate crystals in the pipeline. A considerable amount of hydrate crystal can cause hydrate plugs in the pipeline. The hydrate plug effect lengthens production outages and results in loss of money for the maintainer [1], because elimination of the plug is a time consuming procedure. One of the most widely used traditional solutions to prevent hydrate formation is addition of methanol to the steaming gas. The methanol helps to dehydrate the gas, thus, the growth of hydrate crystal is limited [2]. The methanol can be transported from the gas separator station or gas collection station to the gas wells via a dedicated pipeline. This traditional solution was the most popular practice for a long time in gas industry. This technology is safe, but has several drawbacks, such as the cost of additional pipe to the gas wells, the cost of the methanol regeneration [3], and methanol contamination of the environment. Regardless the methanol technologies are still major solutions for hydrate inhibition for the currently operating gas wells. The modern installations use methanol in low concentrations for newly installed gas wells. These bring the need for injecting methanol locally at the site of the gas wells. Thus, an injection unit is needed for this purpose. The production-related aspects such as the consumers expecting increased flexibility from gas provider cannot be neglected, because of the habits of the consumers and the appearance of competitive energy-saving technologies. Therefore the different injection systems spread noticeably dynamically in recent years. In the first part of the paper, the main features of the commercial injection systems are detailed followed by the introduction of a newly developed injection system. Chemical injection systems In gas industry there are two applications, where chemical injection is required. The first one is the above described gas hydrate inhibition technology. The second one is where usage of corrosion prevention chemicals is needed. New chemical injection technology appeared in sophisticated form using a wide service range as a result of a change in gas production and distribution in the last decade. Structure of pneumatic chemical injection pump systems The most commonly used chemical injection units in gas industry employ the energy of gas for dosing. These systems are called “operating without auxiliary power systems”. In details, the energy of gas steam in the pipeline is used as a power source in the dosing of hydrate inhibiting chemical pump systems. This method has the major advantage that the gas is always available, or it is not necessary to provide a separate power source for the equipment to operate. Big disadvantage of the technique is that gas is required for operation and this gas is emitted to environment, which is a harmful effect, moreover, the cost of it is also not negligible. There are some solutions to utilize the emitted gases, but these methods raise the price of the pump. 8 Fig.1 shows a common structure of chemical injection pump systems, which illustrates how the actuated gas gets into the injection system through a gas preparation unit. This unit is typically a special tank, which is filled with methanol. By using this type of conditioning one can be sure that the gas is sufficiently dry, thus, the water does not disturb the operation of the pump. It is needed to avoid the risk of freezing out various components. The gas bubbles through the methanol, so it loses a significant part its water content. The replacement of the methanol is needed at appropriate intervals of time as a technology management. After the gas leaves the conditioning unit, the gas steams through a two-stage pressure regulator. This unit provides the required pressure to the system for operation. The timer unit is responsible for providing appropriate technological demands of the schedule, i.e. sufficient number of strokes is needed to operate the pump. The chemical is injected with the piston of the injection pump. The injection system includes other components, which are essential for the operation, such as check valves, ball valves and other fittings. The hydrate inhibitor compound can be found in the tank. The main task of the control system is the monitoring of the injection system, so several parameters are required to be measured, e.g. the level of chemical in the tank, the pressure of the after the regulator and the ambient temperature. One of the most substantial information about the injection pump is the entry of the inhibitor to the pipeline. To measure these parameters, transducers are needed, which require some solar energy for the operation. Fig.2 shows a well- organized, modern gas well in Hungary. The injection system used employs the gas of the well for its operation. Injection system from HASKEL MILTON ROY The HASKEL MILTON ROY company develops and manufactures various injection systems for special oil and gas applications for many years. Chemical Injection Pump (CIP) series is one of the most widely used equipment for gas wells. A pneumatic version of the device is shown in Fig.3, as a typical system. The driving gas or air is fed from the regulator to the gas- drying unit. The air is let into the injection pump, which includes the timer unit. The chemical enters the pump through the appropriate valves and filters. The dosage of chemicals can be adjusted through the rate setting gauge unit (typically it is a burette) at, which the quantity of chemical of one injection is displayed. The strokes can be counted by using the flow meter on the system. The check valve prevents the injection pump from the backpressure of the process fluid. The company can deliver the complete unit in different versions according to the quantity of injection, e.g. the given series can be adjusted from 0.27–25.8 dm3 min-1. The range of operation pressure is also wide, typically 120–600 bar. Structure of electric injection systems This category includes systems, where the electricity is the only energy source. In general, the site of the gas Figure 1: Block scheme of chemical injection system Figure 2: Chemical injection system on a gas well Figure 3: Pneumatic injection system from HASKEL MILTON ROY [4] 9 well is not equipped with electric power because of the relatively large distance from the processor (separator station, collector, etc.). Thus, these systems are not commonly employed in the gas well applications. Most of the electric chemical injection pump systems are based on the principle of displacement pumps. These high-pressure pumps provide the chemical injection by using speed adjustment in a wide range, or with different piston volumes. As shown in Fig.4, the electrical injection pump system contains relatively few units. It cannot be ignored that the electric pump does not only include a motor-driven piston with frequency inverter, the controller electronics is also an important part of the system. Solar energy can be applied as energy source of this type of system. This source can feed the control electronic besides the main actuator of the system. In-house designed injection system The development of the injection system was a part of a bigger project with the Scada Ltd, a company from automation industry. The aim of development was to create an injection system, which can be used mainly at Hungarian gas wells. Thus, the temperature requirement of the system was in the -40 °C to 60 °C range. The system must be capable of working in ‘Ex’ environment with high efficiency. The power source of the actuator is solar energy to reach close to zero emission of the system. In the system, corrosive chemical has to be injected, thus, a corrosive resistant is required at critical parts of the equipment e.g. piston, piston space, etc. The maximal working pressure was 160 bars. A wide range of injection volume was also a fundamental requirement. These requirements were the main parameters at the design phase of the equipment. The block scheme of the developed system can be seen in Fig.5. The base of the system is an Ex class asynchronous motor (Table 1), which is assembled with a helical gear unit from Nord manufacturer. The base of injection system is a crank mechanism, which transforms the rotary motion to alternating motion. The pump has two pistons to reach higher injection volume and a well- balanced load of the motor. The main parameters of the mechanics can be found in Table 1. The sealing of the mechanism is a critical design point. The relatively high pressure, corrosive fluids, and high stroke number require special solution for sealing. A single acting, spring energized-seal was selected. The application conditions demanded the tight guiding of neck of piston, to avoid the damage of seals during the operation of pump (Fig.6). The electromechanical actuator gives the opportunity to smoothly alter the injected volume in continuous or in periodic operation. Table 2 shows the injected volume of the pump at different frequency of the inverter in continuous mode. The power source of the system is an array of solar cells. The power of solar energy system was 2 kW to get extended service and short recharge time. The capacity of the accumulators was 800 Ah. The developed equipment has a PLC-based control system to monitor the main parameters of the technology and to change the behaviour of the injection system. Figure 4: Block scheme of electric power injection system Figure 5: Block scheme of new injection system Table 1: Main parameters of the injection pump designed Asynchronous motor Nominal voltage 230 V(Y) /400 V (D), 50 Hz Nominal power 0,12 kW Ex class PTB 07 ATEX 1058 X II2G Ex de IIC T4 Type DEx 63K/4K Vendor HEW Helical gear unit Gear ratio 41,58 Ex class Ex II 2G c IIC T4 X Type SK02F-IEC63 /26-63 S/4 TF/2G Vendor Nord Crank mechanism Length of stroke 20 mm Ratio of mechanics 7 Injection unit Diameter of piston 5 mm Vendor of seals Trelleborg Type of seals Turcon® Variseal® 10 The measured parameters are the pressure of siphon of well, drill pipe, injection pipe; the temperature of pipe of well, soil, chemicals, controller, temperature of inverter; the level of the chemical in the tank; the current and the voltage of the inverter. The frequency of the inverter can be changed in function of the conditions of the well. The injection system is able to communicate wirelessly with a server. By using modern computer technology, the main parameters of the system can be monitored and changed. Fig.7 shows the website of the system, which shows, on one hand, the main actual online parameters of the system, on other hand, the trends and history of the parameters on diagrams. The structure of the website is user friendly, simple. A model technology for an actual gas well was set up at the site of the Scada Ltd. near to Hajdúszoboszló (Fig.8). At the test site, main parameters like level of the chemical in the tank, process gas pressure, frequency of the inverter, etc. were monitored and measured. After the successful test of the model technology, the equipment was moved to a real gas well at Békásmegyer in October 2012. Conclusion A new solar energy-based chemical injection pump system was developed with a modern control and computerized support system at University of Miskolc, Research Institute of Applied Earth Sciences. The equipment meets the modern requirements. It is an efficient injection system; it has wide operation range and has user-friendly high-level computer support. In the process automation field of gas industry, the adaptation of new approaches appear slower in consumer or in industrial environments than in other fields of industry. The new, energy efficient, solar supplied, standalone chemical injection systems hopefully can enter in the field and will be utilized more wells in Hungary. Acknowledgement This work was partially supported by the framework of the Centre of Excellence of Mechatronics and Logistics at the University of Miskolc. REFERENCES [1] BERECZ E., BALLA-ACHS M.: Gázhidrátok, Akadémiai Kiadó, Budapest, Hungary, 1980 (in Hungarian) [2] HAMMERSCHMIDT E.G.: Formation of gas hydrates in natural gas transmission lines, Ind. Eng. Chem. 1934, 26(8), 851 [3] LEDERHOS J.P., LONG J.P., SUM A., CHRISTIANSEN R.L., SLOAN, E.D. JR.: Effective Kinetic Inhibitors for Natural Gas Hydrates, Chemical Engineering Science, 1996, 51(8), 1221–1229 [4] Williams Chemical Metering Pumps: www.willamspump.com (last accessed: May 25, 2014) Figure 7: Website of pump system Figure 6: Section of injection unit Table 2: Injected volume of pump at 150 bars Frequency of inverter, Hz Injected volume, dm3 h-1 30 21.71 50 36.19 70 50.67 90 65.14 Figure 8: Block scheme of chemical injection system HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 11–15 (2013) IONIC LIQUIDS IN ADVANCED ENERGY STORAGE CELLS ATTILA GÖLLEI ! AND ATTILA MAGYAR Department of Electrical Engineering and Information Technology, University of Pannonia, Egyetem Str. 10, Veszprém, H-8200, HUNGARY !E-mail:gollei.attila@virt.uni-pannon.hu The interest in the research and application of ionic liquids has been increased worldwide, especially in petrochemical industry, nuclear industry, oil and gas industry, chemical, and electro-chemical industry. In the field of chemical operations in microwave field with ionic liquids (organic chemical synthesis, catalytic operations, etc.) new horizons opened that were supported by the dramatically increasing number of studies dealing with this question and their impact. Due to the high durability of ionic liquids as a function of temperature, a wider temperature range of operations can be accessed. In addition, it offers environmental friendly solution in the replacement of the toxic solvents with generally low volatility. The area of application extends to electrolytes of novel battery cells. The important prerequisite for application of ionic liquids is the knowledge of their physical and chemical properties. Our aim was to test the dielectric properties, viscosity, and temperature dependence of the electrical conductivity. Considering our measurements, it can be declared that the aforementioned properties show significant temperature dependence in the case of ionic liquids. This knowledge is important for the usability, design and execution of production and related optimization processes. It is particularly important that energy storage cells can be exposed to large temperature fluctuations. This study discusses the sample materials, the usage possibilities of the sample materials and the obtained results. Keywords: Ionic liquids, energy storage cells, conductivity, temperature dependence, viscosity Introduction The key parameters in the modelling of microwave behaviour of ionic liquids are the dielectric constant, dielectric dissipation factors, in addition the parameters characterizing the polarizability and the microwave energy absorption. Traditional measurement methods of the static dielectric constant of ionic liquids fail, because they cause intense shortcut due to their high electrical conductivity. According to literature, these contexts have not been studied systematically because of the novelty of this subject. The significance of these results is critically important in the planning and controlling chemical reactions. Ionic liquids (IL) opened new technological possibilities at the end of the 1990s providing a wide range of applications, especially for areas of chemistry and particularly in green chemistry. ILs are beneficial due to their special chemical and physical characteristics. Nowadays, the use of energy storage units, particularly mobile phones, electric vehicles, and uninterruptible power supply systems has more and more importance. Recently, room temperature ionic liquids (RTILs) have been extensively studied as electrolytes of lithium ion batteries from cellular phones to electric vehicles. The unique properties make RTILs attractive candidates for the electrolyte base of safe lithium battery [1, 2]. Generally, RTILs are well known to be thermally stable and non-flammable and their nature can hopefully improve the safety of electrochemical devices with aprotic solvents, such as Li batteries and super- capacitors. 1-Ethyl-3-methyl-imidazolium tetrafluoroborate (EMIM-BF4) has a high ionic conductivity, which is comparable to those of organic solvent electrolytes. In addition, the viscosity of EMIM-BF4 is low. Therefore, EMIM-BF4 is expected to be a good electrolyte candidate for Li batteries. We have found that a Li/LiCoO2 cell with RTILs as an electrolyte base works reversibly, indicating that QAimide RTILs are quite stable even at the Li reduction potential. Clearly, RTILs improve the safety of Li-ion batteries with carbon anodes. To date, RTILs are the most promising electrolytes for improving the safety of Li-metal batteries, which could be important for higher energy densities [3]. A set of four imidazolium ionic liquids (solid at room temperature) and one imidazolium ionic solid were screened for their potentials as electrolytes in manganese dioxide free LECLANCHÉ batteries, equipped with a zinc anode and graphite cathode [4]. It should be noted that the ionic liquid 1-ethyl-3- methylimidazolium hydrogensulphate (EMIM-HSO4) similar to many other ionic liquids, such as those based on nitrate and dihydrogenphosphate anions, have also been found they work well in this battery design. Novel 12 batteries are designed using standard cathode materials such as MnO2, PbO2, NiO and AgO, and anode materials such as Zn, Sn and Pb. Additionally, by using a solid polymer electrolyte composed of polyvinyl alcohol and anionic liquid, new types of solid state batteries were demonstrated with discharge voltages ranging up to 1.8 V, depending upon the type of cathode and anode used [5]. Ionic liquids like 1-butyl-3-methylimidazolium tetrafluoroborate (IMIM-BF4) or hexafluorophosphate (IMIM-PF6) and 1-butyl-4-methylpyridinium tetrafluoroborate (PyBF4) were mixed with organic solvents such as butyrolactone (BL) and acetonitrile (ACN). The specific ionic liquid used in the given study is shown in Fig.1. A lithium salt (LiBF4 or LiPF6) was added to these mixtures for possible application in the field of energy storage (batteries or super-capacitors). Viscosities, conductivities and electrochemical windows at a Pt electrode of these electrolytes were investigated. All the studied electrolytes are stable toward oxidation and exhibit a vitreous phase transition, which has been determined by application of the conductivity measurements. Mixtures containing the BF4 - anion exhibit the lowest viscosity and highest conductivity [6]. Methods The essential precondition for the application of ionic liquids in batteries is a preliminary examination of its properties. Therefore, every measurement method available to us was used on a number of ionic liquids examined. We have carried out microwave measurements, such as temperature change, dielectric properties at 2.45 GHz frequency, electrical conductivity according to temperature and viscosity changes depending on the temperature. The velocity of temperature change of ionic liquids was measured in the CEM Discover unit. This is a well- known and commercially available apparatus in microwave chemistry. The device has a cylindrical operating space and on the cylindrical peripheral surface there are many slots where microwave energy can enter. This solution provides a special high homogeneity of the microwave field. The temperature was measured at the bottom of the compartment with an infrared thermometer. The amount of the samples used was 0.5 g, which was placed in an inner diameter of 12.5 mm of borosilicate cylindrical glass flask. It is important to know that during microwave treatment the microwave energy interacts with the treated material, and the material converts the electric energy into thermal energy according to its characteristic of dielectric properties. During this conversion process, the macroscopically observable result is the measurable increase of the temperature in the treated material. The rate of the temperature increase depends on the microwave field and the treated material properties, which is described in the following equation: ΔT /Δt = Pv / ρCp = iE 2 fε"/ ρCp (1) where, the ∆T/∆t is the velocity of temperature change in Kelvin per second, Pv is the absorbed power in the sample, ρ is the density of the sample, Cp is the specific heat of the sample, E2 is the strength of the electromagnetic field in the sample, f is the frequency of the field and ε” is the dielectric loss, i means this is a complex value [4]. Eq.(1) shows that the temperature of the treated material is influenced by a number of factors. The separate testing of their effects is not an easy task, since ρ, ε, Cp features themselves are temperature dependent, and this value is difficult to accurately measure inside a material. In case of strictly homogeneous series of examined compounds the situation is more simplified, because some simplifications are allowed. Assuming ρ and Cp does not change significantly in the function of temperature, the multiplication ρ Cp was almost considered to be constant, so that the rate of temperature change is determined by E and ε" only. Further simplifications can be made if the device is single-mode with the same volume and shape, and the microwave energy is constant during the investigation. The value of E is determined by ε' and ε' and ε" and the rate of temperature change is essentially determined by ε" value. Figure 1: The structure of investigated ionic liquids Table 1: The formulas of investigated ionic liquids Compound R Y 1 CH3 BF4 2 C2H5 BF4 3 C3H7 BF4 4 C4H9 BF4 5 CH3 PF6 6 C2H5 PF6 7 C3H7 PF6 8 C4H9 PF6 9 C2H5 CI 10 C2H5 Br 11 C2H5 SCN 12 C2H5 N(CN)2 13 C2H5 N(SO2CF3)2 13 Measurement setup The schematic representation of the in-house designed experimental set-up for automated online measurement of dielectric properties of ionic liquids in a definite temperature range is shown in Fig.2. It is composed of the following devices and instruments; cylindrical sample holder unit, thermostat, peristaltic pump, waveguide, temperature sensor, displaceable piston, stepper motor, magnetron, detectors, control unit, and a PC. The IL sample, which is placed in the thermostat to be kept at the desired temperature, is passed across the waveguide having a length of about 3λ through the sample holder tube with the help of a peristaltic pump. The electric energy is transformed into microwave energy by the magnetron. The stepper motor is controlled by a microprocessor control unit equipped with an Intel 8-bit microcontroller, a 12 bit A/D converters for receiving the four diode-detector signals, a stepping motor driver, and a RS-232 serial interface to connect it to a PC. The control unit collects the detector signals and the temperature data from the temperature sensor and sends them to the PC. Furthermore, it controls the position of the short circuit displaceable piston on the basis of algorithm software elaborated for this purpose. The measurement method is based on the compensation of phase change due to the microwave energy absorption of the liquid sample. The short circuit piston situated behind the sample must be actuated for compensation. The energy conditions created by the wave front in the waveguide are measured by four diode-detectors. Dielectric constants, dielectric loss factors and the temperature dependence of the dielectric properties of ionic liquids intended to be used in batteries were determined by the above described self-designed microwave dielectrometric apparatus (Fig.3) at the frequency of 2.45 GHz and at different temperatures (30°C, 40°C, 50°C, 60°C, 70°C, 80°C, 90°C, 100°C, 110°C and 120°C). The speed of the change in temperature depends on the electrical field strength in the material (E), the absorbed microwave power, (Pv) density (ρ), the specific heat capacity (Cp) and the dielectric loss factor ε” and can be given by Eq.(1) [7]. Results and Discussions The effects of R groups (Table 1) were examined in case of two anions - BF4 (1–4) and PF6 (5–8) at 3 to 5 W of microwave energy. The results indicate that the decrease of the number of C atoms in R increases the rate of temperature change. By increasing the microwave power significantly the temperature rapidly converge and sometimes crossing over, but this is not typical in Fig.3. Overall, the rate of temperature rising follows the order of ε' and ε" suggesting that in a strictly homologous series the temperature change of ionic liquids is determined by ε" [5]. After evaluating the rate of temperature change, we investigated two additional properties of the ionic liquids used in the advanced batteries. These parameters were the viscosity and the electrical conductivity. Both of these parameters are the functions of temperature. They were examined between 20 °C and 100–120 °C temperature range. The viscosity was measured up to 95 °C due to the sample material being heated by a water- filled heat exchanger. The measurements were performed in a SV-10 type of vibration viscometer (A&D Ltd. Japan). Fig.4 shows that the viscosity at room temperature is high, typically several hundred Pa·s. It decreases rapidly with the temperature increase similar to a reciprocal curve. At 100 °C, the value of viscosity can be ten times lower then at the initial value as can shown in Fig.4. The investigation of the viscosity is important, since it basically correlates with the mobility of ions in the electrolyte, which can affect the conductivity, and thus the internal resistance of the battery cell. By cooling ionic liquids, the viscosity is increasing rapidly. Around freezing temperature the ionic liquids were already similar to a “honey density” mass. Since the energy storage cells are used in a wide temperature range, it is important to know that the Figure 2: Scheme of the experimental setup Figure 3: Temperature profile for 1-butyl-3- methylimidazolium tetrafluoroborate at 3W 14 electrical conductivity of the electrolytes is a function of temperature. Fig.5 contains five different electrical conductivity plots of ionic liquids in the function of temperature. It can be seen that the above mentioned BMIM-BF4 does not have the highest conductivity value, but the depending on the temperature the conductivity value of this material increases with the greatest intensity [6]. A new and surprising result can be seen in Fig.6. The dielectric constant value of BMIM-BF4 reaches a value of 1 at close to 100 °C, similarly to the value of the vacuum and with the further raising of the temperature this value decreases into the negative range. This fact suggests that by increasing the temperature the electrical conductivity is rising greatly and this compound is no longer an insulating material, but a conductor. Therefore it is important to keep the temperature in an adequate range when batteries are used and especially in the charging period. The dielectric constant (E1), dielectric loss factor values of ILs built up of the same BF4 - anion and six previously described cations containing alkyl chains with different lengths (DiEtMeIm+, DiEtEtIm+, DiEtPrIm+, and DiEtBuIm+) at different temperatures between 30°C and 120°C are shown in Fig.7. At the initial measuring temperature of 30 °C all four of the studied ILs have similar dielectric constants around 7. With the increase in temperature the E1 values for DiEtMeIm(BF4), DiEtEtIm(BF4) and DiEtPrIm(BF4) slightly increase up to 13, while the dielectric constant for DiEtBuIm(BF4) shows a sudden break at 90°C and at 120°C it reaches the value as high as 27. This could be explained with some sudden changes in the structure of the ILs or in the physico-chemical interactions between the anion and the cation. Excluding the results for DiEtBuIm(BF4) the values at elevated temperature show that the highest E1 value belongs to DiEtMeIm(BF4), followed by DiEtEtIm(BF4) and then by the DiEtPrIm(BF4), hence the dielectric constant increases with the decrease in the alkyl chain length of the cation. Taking into account the dielectric constant, dielectric loss and electrical conductivity results of all the investigated ILs it can be concluded that the alkyl chain length of the cation and the structure of the anion strongly influence the dielectric properties of the ILs and that the highest G value is exhibited by the DiEtMeIm(BF4) IL at 30°C temperature, hence it is the most suitable candidate for battery applications. Furthermore it can be stated that implementing the knowledge about the connection between the IL structure and the dielectric properties another ILs should be studied in order to select the ones most adequate for electrolyte application. Figure 4: Viscosity of BMIM-BF4 as a function of the temperature Figure 5: Conductivity of five compounds as a function of temperature Figure 6: Dielectric constant of BMIM-BF4 as a function of the temperature Figure 7: Dielectric constant of ILs containing common BF4 - anion and different cations 15 Conclusion Due to the special characteristics of ILs, such as wide electrochemical windows, high inherent conductivities, high thermal and electrochemical stability, tuneable physicochemical properties, etc., they are potentially excellent candidates for environmentally sound, green electrolytes in batteries. In order to predict their success in a specific application, it is essential to gain information about their dielectric properties. However, several techniques have already been published about determination of the dielectric properties of ILs, so far no experimental methods were introduced due to the fact that such measurements cannot be realized by classical methods because of the high conductivity of ILs. Summarizing our results, we can conclude that when ionic liquids in the energy storage cells were used, their physical, chemical and electrical parameters are significantly temperature dependent. In any case, when they are applied, we should specify the range of application temperature. Stepping out of this range of the energy storage cell may not meet the expected specification values, or in the worst case it may be permanently damaged. Acknowledgements We acknowledge the financial support of the Hungarian State and the European Union under the TAMOP- 4.2.2.A-11/1/ KONV-2012-0072. REFERENCES [1] ALARCO P.J., YASER A.L., RAVET N., ARMAND M.: Lithium conducting pyrazoliumimides plastic crystals: a new solid state electrolyte matrix. Solid State Ionics, 2004, 172(1–4), 53–56 [2] SCHIFFMANN R.F.: Microwave and dielectric drying, in: 2nd Handbook of Industrial Drying, (Ed. MUJAMDER A.S.), Vol 1, Marcel Dekker, New York, USA, 1995, pp. 345–372. [3] HIKARI S., HAJIME M., KUNIAKI T.: Application of room temperature ionic liquids to Li batteries Electrochimica Acta, 2007, 53, 1048–1054 [4] ZHANG Z., GAO X., YANG L.: Electrochemical properties of room temperature ionic liquids incorporating BF4 and TFSI anions as green electrolytes, Chinese Science Bulletin 2005, 50(18), 2009 [5] GAO X.H., ZHU M.J., ZHANG Z.X.: Research on room temperature molten salts with alkylated imidazolium salt as electrolytes, Chemical World (Supplement, I), 2004, 45, 156–157 (in Chinese) [6] NISHIDA T., TASHIRO Y., YAMAMOTO M.: Physical and electrochemical properties of 1-alkyl-3- methylimidazolium tetrafluoroborate for electrolyte, J. Fluor. Chem., 2003, 120, 135–141 [7] GÖLLEI A., VASS A., PALLAI E., GERZSON M., LUDÁNYI L., MINK J.: Apparatus and method to measure dielectric properties (ε’ and ε”) of ionic liquids. Review of Scientific Instruments, 2009, 80(4), 6734–6748 Advertise upcoming meetings, conferences, and workshops; make public announcements; introduce your research laboratory; a new product or a service in the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-298, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 17–26 (2013) MONITORING AND DIAGNOSIS OF MANUFACTURING SYSTEMS USING TIMED COLOURED PETRI NETS ADRIEN LEITOLD1, BRIGITTA MÁRCZI2, ANNA IBOLYA PÓZNA2, AND MIKLÓS GERZSON2! 1 Department of Mathematics, University of Pannonia, Egyetem str. 10, HUNGARY 2 Department of Electrical Engineering and Information Systems, University of Pannonia, Egyetem str. 10, HUNGARY !E-mail: gerzson@almos.uni-pannon.hu Novel fault modelling and integration method were applied in the case when the faultless operation of the system was modelled by a high-level, coloured Petri net. In order to achieve realistic investigations, a timed coloured Petri net model of the system was constructed, where faults can occur in the manufacturing lines. The faultless and fault containing models were implemented in CPNTools both for non-timed and timed cases. The resulted model was investigated both via simulation and using the occurrence graph. For efficient analysis of the occurrence graph a software module called OGAnalyser was developed. Keywords: monitoring and diagnostics of manufacturing processes, coloured timed Petri nets, probabilistic models, occurrence graph Introduction Models are often used for the description and investigation of complex systems even if they cannot perfectly describe the investigated system. The course of a manufacturing system can be split up to distinct steps of serial or parallel technological sub-processes. This enables the description of the system by a discrete event systems model [1] in the form of Petri nets. During the design process of a manufacturing system, often only a model of the faultless operation is created. The integration of possible faults into the model could give important information for more complex investigation of the system. In our previous work [2] we have integrated fault events with different occurring possibilities into low-level Petri net models of manufacturing systems in such a way that the size of the model remained almost the same. In our recent work, we applied the above fault integration method for the case when the faultless operation of the system was modelled by a high-level, coloured Petri net (abbreviated as CP-nets) [3]. In low- level Petri nets the transitions fire instantaneously, but the events of a real system take place for a certain amount of time influencing the operation of the system. Therefore, a timed coloured Petri net model of the system is constructed in order to achieve realistic investigations. CPNTools [4] offer tools for modelling and analysing of CP-nets. There are two possibilities for the investigation of a manufacturing system in CPNTools: the simulation and the analysis of the occurrence graph. In case of fault modelling using different occurring possibilities; however, the standard occurrence graph does not give information about the probability of the different occurring states of the system. Therefore, a special software module, called the OGAnalyzer has been developed for solving this problem. Weights can be assigned to the arcs of the occurrence graph and the software calculates the probability of each node in the occurrence graph, i.e. of each operational state of the system. Petri net model of a manufacturing process Petri nets enable both the mathematical and the graph representation of a discrete event system to be modelled, where the signals of the system have discrete range space and time is also discrete [5]. Petri nets can be used for describing a controlled or open loop system for modelling the events occurring in it, and for analysing the resulted model. For the different application purposes, various modifications of the original Petri net were developed with the aim of improving the modelling capabilities. One of the approaches is the coloured Petri nets (CP-nets). We use here the CP-nets for modelling technological systems and their diagnosis, i.e. for the determination of faulty operational modes of the investigated system. Coloured Petri nets The CP-nets combine the modelling advantages of Petri nets and the compactness of the functional 18 programming language Standard ML [6]. A Petri net is bipartite graph having circles and rectangles as nodes. Circles refer to the ‘places’ in the net and rectangles to the ‘transitions’. Places represent the state of the elements in the modelled system, while transitions correspond to the actions taken place in it. There are ‘arcs’ between places and transitions referring to logical relations of the system. If an arc directs from a place to a transition then the place acts as a precondition of the given transition while the arcs in the opposite direction represent consequences of transitions. Each place can be marked with one or more coloured tokens representing the state of the modelled element. Here we emphasize two important novelties of CP- nets only. The tokens describing the state of the system have data value, the so-called token colour attached to them. In this paper these colours are used to identify workpieces and to describe the operation to be performed on them. Places, transitions and arcs can have ‘inscriptions’. An inscription of a place determines the set of colours that a token on the place can have. Another place inscription gives the actual number of the tokens on the place, i.e. the current marking of that place. The inscriptions of transitions can contain different types of functions. These functions determine the type of the colour set of the incoming and outgoing tokens and the operation performed on them. The arc inscriptions can be used for evaluating the result of the performed action at the previous transition. These conditional expressions define the colour of the token on the following place. The original Petri net concept did not contain the time; however, the firing of transition takes place instantaneously. In case of real technological systems the time has a great role during the occurring of events. CP-nets also offer the possibility of adding time to the operation of transitions. The firing rule of a transition in timed Petri nets is as follows: a transition is enabled if all of its input places contain sufficient number coloured tokens defined by the arc functions and the allotted time has elapsed. One of the main advantages of modelling with Petri nets is the ability of describing sequences of discrete events that occur both in a serial and in a parallel way. In case of parallelism, we can distinguish two different situations. In the first case the two or more series of events can take place independently of each other. This situation occurs when workpieces can be elaborated in different manufacturing lines in parallel way. In the other case only one of the event sequences can take place because these events exclude mutually each other. These events have the same precondition, and the occurrence of any of them makes this precondition invalid. This kind of parallelism is called ‘conflict situation’. In a Petri net the conflict can be recognized when a place is the precondition of two or more transitions. In this case it is randomly selected, which transition takes place. The conflict can have two different sources in technological systems. A conflict occurs in a technological sense when two or more processes want to use the same tool or resource e.g. a robot. Usually it is worth to assign priority to each conflicting transition in order to define their sequence. When a fault occurs during the operation of the system, it also causes a conflict situation. This can be avoided by adding a special probability function to arc expression functions. By evaluating this function, the occurrence of the fault can be unambiguously determined during execution. Analysis of Petri nets There are two basic directions for the analysis of a Petri net: (i) the structural method, which is independent of the initial state of the net and (ii) the investigations based on a given initial state (the behavioural analysis). In this paper the latter is used for the investigation of technological system behaviour. Simulation is our primary tool for the checking the correctness of a model. Starting from a given initial state the user can check whether the operation of the system terminates in the appropriate state. It can also be investigated, which transitions become enabled in certain steps, whether there is a conflict among then. Simulation investigations do not give unambiguous answers to questions of formal analysis formulating in Petri net literature [7] but they complement them well. Another Petri net analysis method uses the ‘occurrence graph’. The basic idea of the occurrence graph is to construct a graph, which contains all of the reachable markings from a given initial state. These marking are the nodes of the occurrence graph and the arcs connecting the nodes refer to the logical relations realized by the firing transition between two markings. Unfortunately, the occurrence graph even of a small Petri net may become very large. Therefore, several reduction methods were proposed in order to get a relatively small occurrence graph [7]. Most of the simulation tools, as the CPNTools [4] used by us, are able to construct the occurrence graph. Modelling and analysis of technological systems In the following, we wish to demonstrate the use of CP- nets and their occurrence graphs for modelling and analysis of technological processes. Both for normal faultless and faulty mode operation of the technological system are considered in the non-timed and timed cases. The analysis is performed via simulation and with different investigations of the occurrence graph. The manufacturing system and its operating procedure A simple case study is presented here for a manufacturing system containing two manufacturing lines and a robot. The arrangement of the system can be seen in Fig.1. 19 The workpieces to be processed appear on the input place IN. The task of the robot is to put them to the appropriate input place of a manufacturing line M1_IN or M2_IN according to operational instructions. Assume that the two manufacturing lines perform different actions on the workpiece of interest. When the manufacturing process is over, the finished workpiece appears at the output end of the line either on M1_OUT or M2_OUT depending on the performed action. If the workpiece has to be modified on the other manufacturing line then the robot puts it onto the other input place. If the manufacturing process is over then the robot puts it into the product container OUT. Assume that one workpiece at a time can be on the input and output places of manufacturing lines. It follows that the robot can only transfer a workpiece from place IN to the input place of a manufacturing line if this place is empty and the precondition of the start of a manufacturing process is that the output place of this line should be empty. As a general precondition of all transfer processes the robot has to be free. The CP-net model of the normal (faultless) operation of the manufacturing system in the form of a screenshot from CPNTools can be seen in Fig.2. The process starts with a token at place START. The transition generator generates the prescribed number of tokens representing workpieces at place IN. The colour assigned to a token contains an identifier of the workpiece and a code referring to the manufacturing process or processes to be carried out. Four kinds of manufacturing mode are possible in this manufacturing system: the workpiece has to be processed on line 1 (denoted by m1) or on line 2 (m2) only, or it has to go through the line 1 and then through line 2 (m12) or in reverse order (m21). As an example, the token (1,m12) refers to the workpiece having identifier '1' and this piece has to be processed first on line 1 then on line 2. The state of input and output of manufacturing lines is modelled with two places. The places Tin_empty, Tin_full refer to state of inputs and they are mutually exclusive. The colour of tokens referring to the state of these places consists of the identifier of the line only. If the input place of the line 1 is empty then there is a token having colour m1 on the place Tin_empty, and there has not to be a token having colour m1 on the place Tin_full. The same applies for the places (Tout_empty and Tout_full) describing the state of the output places of manufacturing lines. Figure 2: The CP-net model of the faultless operation of the manufacturing system Figure 1: The manufacturing system 20 The transfer processes are disintegrated into three steps: to transitions referring to the (i) start and (ii) end of the transfer and to a place representing the (ii) transfer itself. Transitions Transfer#_start correspond to the start of transfer processes from IN, if # = 1, from Tout_full if # = 2 or 3. The places Transfer# represent the transfer process 1, 2 or 3. The transitions Transfer#_end refer to completing of transfer process to the input place of a line (# = 1), to product container (# = 2) or to the input place to another manufacturing line (# = 3). The place Manuf_line refers to two manufacturing lines and the colour of token shows the line being processed. Assume that only a single fault can occur in the system during manufacturing: the identification label of the piece can get damaged therefore it cannot be identified. Workpieces with damaged label get into a separate container represented by place OUT_Fault. The repairing of the label is not handled in this example. The modified part of the Petri net model can be seen in Fig.3, where the occurring of fault is taken into account. The occurrence of the fault is forced by a check function built into arc inscriptions in the Petri net model randomly. This check function returns with a fault in predefined probability. This probability of the can be set in the definition part of the net and different fault probability values can be assigned to the two manufacturing lines. By comparing Figs.2 and 3, it can be stated that a new place OUT_fault appears as a consequence of fault modelling and integration and the arc expression functions of arcs starting from transition Manuf_line are extended with the fault checking. For realistic investigation of a technological system the timed version of the CP-net model is used as a case study. A time point is assigned to the transitions. The transition generator does not belong to the technological system closely, so it fires under zero time, i.e. instantaneously. Different time units are assigned to the other transitions. These time values appear as transition inscriptions ‘@+i’$ (where$ i is an integer number defining the amount of time in seconds) on the net as it can be seen in the Fig. 3. During the simulation investigations different time values have been applied in order to check the possibility of a deadlock. Preliminary analysis by simulation As a first step, simulation is applied for the investigation of the developed Petri net model of the manufacturing system. These were carried out assuming both faultless system operation and when fault can occur during the manufacturing. Both non-timed and timed operational modes were considered. The short description of the simulation is as follows. The simulation starts with the generation of tokens representing workpieces. The number of these tokens i.e. the number of workpieces to be processed can be modified in the arc expression function belonging to the transition generator. The type of processes to be performed i.e. the label referring to the manufacturing line(s) of the tokens can be set in the definition part of the CP-net. As a result, the colour of generated tokens refers to the identifier and to the process(es) to be carried out. Because of the highest priority of transition generator the transition Transfer1_start can only fire after the prescribed number of token appears at the place IN. As mentioned above, the transition Transfer1_end refers to the completion of the transfer of the workpiece to the input place of a manufacturing line. As a next step, the firing of transition Manuf_line corresponds to the completion of manufacturing process. The only fault in the system can occur during this process. If it happens then the token gets into the place OUT_fault, which represents the fault container. If the fault does not occur then the system checks whether the manufacturing process is over or the workpiece has to be also processed on the other line. In the first case the robot puts it to the place OUT representing the product container, while in other case it transfers the token representing the piece into the input place of another line. The modelling of the other two transfer processes (transfer of a workpiece from the output place of a manufacturing line either to the product container or to the input place of the other line) is similar to the transfer from place IN to place Tin_full. If the appropriate input place and the robot are free, another transfer process can start. The primary goal of the simulation is to check the correct operation of the model. Another aim is the investigation of possible deadlock situations. In case of deadlock, the process stops without completing all the prescribed technological actions. It can happen in case of timed simulation of our investigated system and it refers to the wrong determination of the timing of actions. Figure 3: The Petri net model after integration of fault 21 Simulation cannot ensure the thorough investigation of the modelled system, but it complements well the further analysis. For the detailed investigations we applied the analysis of occurrence graph. Analysis based on the occurrence graph The thorough analysis of the behavioural properties of a CP-net can be performed using its occurrence graph. The concept of the occurrence graph was introduced above. The CPNTools generates automatically the occurrence graph, but the check functions used for the fault generation have to be removed from the arc expressions otherwise the occurrence graph is generated only for the normal mode or for the faulty mode. For the illustration of generation and analysis of the occurrence graph let us assume that there are two workpieces to be processed, one of them has to go through manufacturing line 2, while the other has to go first through manufacturing line 1 then through line 2. Let the operation of the system be faultless and let the firing of all transition be instantaneous, i.e. the net is non-timed. The resulted occurrence graph can be seen in Fig.4. The explanation of numbers in the occurrence graph can be seen in Fig.5. The frames in the upper part of boxes are the identifiers of system states. The expression x:y in lower part shows the number of preceding and succeeding states. The token distribution belonging to a node can be obtained by selecting the triangle in left low corner. The opening window contains the name of places and the number and colour of tokens. All the branches on the graph are explained by technological reasons, as there is no built in rule for the robot to start the transfer with any particular workpiece. The graph has only one terminal node (highlighted by green) and it refers to the normal termination of the process i.e. all the prescribed manufacturing processes terminated properly. Now we can repeat the simulation with the same initial condition relating to workpieces but assuming that faults can occur during manufacturing process. The resulted occurrence graph can be seen in Fig.6. Figure 4: The occurrence graph - no fault, no timing Figure 6: The occurrence graph - with fault, no timing Figure 5: Key for the occurrence graph 22 It can be stated that the occurrence graph has become larger and the number of terminal nodes has increased due the effect of the possible faults. The reason of branches on the graph can be either technological (at the nodes 3, 9, 28, 51 and 54) or due to the fault. Only one of the terminal nodes refers to the normal termination of the process (highlighted by green, node 62), the others belong to the faulty cases (highlighted by red, nodes 22, 37, 43, 56 and 57). In case of faulty operation the identification tag of at least one workpiece get damaged during the manufacturing process. Next we repeat the simulation again but adding the timing information to the net and assuming faultless operation. Let the time duration of transfer processes be equal to unit time, and the time of manufacturing processes is equal to 5 time units. The resulted occurrence graph can be seen in Fig.7. The resulted occurrence graph is a tree and its two terminal nodes (nodes 28 and 29) differ only in time stamp of the tokens (the time stamp can be seen only in CPNTools with a detailed label for each node). Comparing this graph with the occurrence graph in Fig.4, it can be stated that the number of parallel branches is less and they have different terminal states. Adding timing information to the model reduces the number of the possible different technological variants. As the fourth case, let the simulation be performed with both timing and faults. The resulted occurrence graph (see Fig.8) is a tree again, the number of nodes is almost the same but the number of terminal nodes has doubled comparing to the occurrence graph in Fig.6. It can be stated based on Figs.4 and 6–8 that the occurrence graph is relatively simple in case of small number of workpieces and the analysis of nodes can be done manually. Thus, it is easily to find the terminal node or nodes referring to normal, faultless termination of the process, and those terminal nodes where the manufacturing of one piece or of both pieces ends with fault. However, the size of the occurrence graph grows Figure 8: The occurrence graph - with fault and timing Figure 9: The node numbers in occurrence graph Figure 7: The occurrence graph - no fault, with timing 23 exponentially if the number of pieces becomes larger, as it can be seen in the left part of Fig.9 assuming faultless or faulty operational mode, and non-timed or timed net. In case of large number of simpler workpieces when the identification is not necessary for each item, the identification tag can be omitted from the colour of the token. Assuming this situation all the simulation investigations (with and without fault, with and without timing) is repeated. The size of the resulted occurrence graphs is less than an order of magnitude simpler compared to the equivalent case with identifier as it can be seen in the right part of Fig.9. Since the workpieces have to be processed on the same manufacturing line, they have the same colour, thus they are indistinguishable. Since the robot selects among workpieces having the same colour, there will be no technological branches on the occurrence graph. This results in a much simpler occurrence graph especially when there are large number of workpieces and the processes to be performed is one or two kind of sorts. As an example, let us investigate the structure of occurrence graphs in case when there is no timing in the system, the token colour contains the identification tag and both faultless and faulty operational mode is assumed. According to data in Table 1, the size of the occurrence graph depends on the number of workpieces, the number of manufacturing procedures to be processed, and the presence of fault. If fault can occur during the manufacturing process both the number of nodes and the arcs increases dramatically. There are several terminal nodes, too, but only one refers to normal termination of manufacturing processes the others belong to different faulty situations. The labelling of workpieces has also a significant effect on the complexity of the occurrence graph. If there are one or more workpieces, which should be processed on both manufacturing lines then the number of nodes and arcs is doubled as it can be seen in the corresponding rows in Table 1. The results of an investigation with timing information can be seen in Table 2. If the number of workpieces is less than four or the workpieces have to be processed on only one manufacturing line, then the size of the occurrence graph depends on the number of workpieces and the presence of fault. On the other hand if there are at least four workpieces and at least two of them have to be processed on both manufacturing lines, but in reverse order then a deadlock situation can occur. A further condition of a deadlock that the manufacturing time should be longer than the transfer time but it is true in general. In case of a deadlock, the process stops because the precondition of transfer processes cannot be fulfilled. There are workpieces on input and output places of both manufacturing lines and therefore no further steps are enabled. If all the workpieces have to be processed on both manufacturing lines and there is no fault then all of terminal nodes refer to a deadlock as it can be seen in the rows marked by a footnote in the Table 2. The identification tag can be omitted from the colour of token in certain cases. If the number of workpieces is small and they have to be processed in different ways then there is a significant change in the occurrence graph. On the other hand in case of large number of workpieces, the structure of occurrence graph becomes much simpler if the identification tag is removed from the colour. Analysis of the occurrence graph using the OGAnalyzer As mentioned above, CPNTools cannot use the information about the probability of faults at the generation of occurrence graph. However, assigning this value to the appropriate edges, the probability of each node of the occurrence graph, i.e. of each system state can be determined. For this purpose, software called Table 1: Comparing the structure of occurrence graphs – no timing No. of work- pieces label (identifier, manufacturing information) faultless operational mode faulty operational mode nodes arcs TNa nodes arcs TNa 1 (1,m1) 7 6 1 8 7 2 1 (1,m12) 10 9 1 12 11 3 2 (1,m1)++(2,m2) 29 32 1 42 50 4 2 (1,m1)++(2,m12) 37 40 1 60 70 6 2 (1,m2)++(2,m12) 39 41 1 62 72 6 2 (1,m12)++(2,m21) 42 44 1 82 95 9 3 (1,m1)++(2,m2)++(3,m12) 145 173 1 306 409 12 3 (1,m2)++(2,m12)++(3,m21) 165 189 1 414 536 18 3 (1,m1)++(2,m1)++(3,m1) 107 123 1 210 273 8 4 (1,m1)++(2,m2)++(3,m12)++(4,m21) 631 772 1 2063 2900 36 4 (1,m1)++(2,m1)++(3,m1)++(4,m1) 340 404 1 949 1300 16 5 (1,m2)++(2,m12)++(3,m21)++(4,m1)++(5,m2) 2208 2797 1 9813 14420 72 5 (1,m21)++(2,m12)++(3,m21)++(4,m12)++(5,m21) 3259 3839 1 23709 32939 243 5 (1,m2)++(2,m1)++(3,m2)++(4,m1)++(5,m2) 1621 2165 1 5252 8051 32 6 (1,m2)++(2,m12)++(3,m21)++(4,m1)++(5,m2)++(6,m12) 8698 10980 1 62751 93323 216 6 (1,m21)++(2,m12)++(3,m21)++(4,m12)++(5,m21)++(6,m12) 12184 14562 1 146992 208602 729 6 (1,m2)++(2,m1)++(3,m2)++(4,m1)++(5,m2)++(6,m1) 5638 7686 1 21630 38886 64 a terminal nodes 24 OGAnalyzer has been developed. We assumed that the occurrence graph belonging to a given initial state of a CP-net model is finite and acyclic. The occurrence graphs of Petri nets modelling manufacturing systems fulfil this assumption in general. Let the probability of faults be known from technological consideration and let the first step of the analysis be the assignment of the arc weights to the edges of occurrence graph as follows. 1. If a node on the occurrence graph has only one outgoing arc, then the next state follows unambiguously, thus the arc weight is equal to 1. 2. If there is more than one outgoing arcs from a given node then it means that different states can follow from it. These states come into existence with different probabilities depending on whether this branch has technological or fault related reason. a. If there is no fault in the system then every branch has a technological reason, because there is no built-in priority rule for the robot to the selection among the workpieces. It means that every selection that is every arc has the same probability so all of the arcs starting form this node have to get the same arc weight, which is equal to the reciprocal value of the number of outgoing arcs. b. The introduction of the fault into the model causes the appearance of another type of branching in the occurrence graph. Let us assume that only one type of fault can occur in a given system state. It results in 2 new system states: one for the normal operation and one for the faulty mode. Let the probability of fault be equal to Pf. Then the weight of arc leading to faulty mode is equal to the probability of the fault while the arc leading to the normal operational mode gets the value 1-Pf. Assigning these arc weights to the edges of the occurrence graph the probability of a given state on the graph can be determined in the following way if the faults occurring one after the other in the system are independent: 1. If the graph is a tree or only one route leads to the given node then the probability of the state representing by this node is equal to the product of arc weights along the route leading from the node representing the initial state to the given node. 2. If there are more than one route leading to the given node from the initial node then the probability value of each route has to be determined with the production of arc weights along the route as in the first step, then to sum these resulted values. 3. If the faults are not independent from each other, then the probability of nodes can be calculated in a similar way but using conditional probability values. As described above, this determination method and the operation of OGAnalyzer has been illustrated using the example of two workpieces to be processed and one of them has to go through manufacturing line 2, while the other has to go first through line 1 then trough line 2. Faults can occur during the manufacturing process and let the fault probability be equal to 0.3 in case of line 1 and 0.1 in case of line 2. As before, we consider the timed case, when the transfer transitions have the same transition time of 1 time unit, while the manufacturing time is equal to 5 for both lines. After the simulation CPNTools generates the occurrence graph, the structure of, which is the input of OGAnalyzer. As a first step OGAnalyzer reads the structure of the generated occurrence graph from the data file generated by CPNTools and visualizes it in its own window as it is shown in Fig.10. The user can get the token distribution belonging to nodes as it is shown at node 1 in Fig.10. The next step is the identification of branches. The software can distinguish between the two different types of branches on the occurrence graph (that are technological and fault caused branches). For this the user has to define the fault colour for the appropriate branches in a separate window. Table 2: Comparing the structure of occurrence graphs – with timing No. of work- pieces label (identifier, manufacturing information) faultless operational mode faulty operational mode nodes arcs TN/Da nodes arcs TN/Da 1 (1,m1) 7 6 1/0 8 7 2/0 1 (1,m12) 10 9 1/0 12 11 3/0 2 (1,m1)++(2,m2) 23 22 2/0 41 40 8/0 2 (1,m1)++(2,m12) 29 28 2/0 59 58 12/0 2 (1,m2)++(2,m12) 29 28 2/0 58 57 12/0 2 (1,m12)++(2,m21) 35 34 2/0 87 86 18/0 3 (1,m1)++(2,m2)++(3,m12) 73 74 4/0 264 263 48/0 3 (1,m2)++(2,m12)++(3,m21) 86 85 4/0 392 391 72/0 3 (1,m1)++(2,m1)++(3,m1) 85 84 6/0 262 261 48/0 4 (1,m1)++(2,m2)++(3,m12)++(4,m21) 173 172 8/2 1473 1472 258/2 4 (1,m1)++(2,m1)++(3,m1)++(4,m1) 365 364 24/0 2121 2120 384/0 5 (1,m2)++(2,m12)++(3,m21)++(4,m1)++(5,m2) 574 573 24/4 9320 9375 1582/8 5b (1,m21)++(2,m12)++(3,m21)++(4,m12)++(5,m21) 154 153 24/24 12076 12075 2196/60 5 (1,m2)++(2,m1)++(3,m2)++(4,m1)++(5,m2) 514 513 24/0 5152 5247 780/0 6 (1,m2)++(2,m12)++(3,m21)++(4,m1)++(5,m2)++(6,m12) 1725 1732 80/40 71817 72336 12300/116 6b (1,m21)++(2,m12)++(3,m21)++(4,m12)++(5,m21)++(6,m12) 325 324 72/72 45541 45540 8280/432 6 (1,m2)++(2,m1)++(3,m2)++(4,m1)++(5,m2)++(6,m1) 2341 2412 72/0 39349 41076 5760/0 a terminal nodes/deadlocks, b labels where all terminal nodes refer to deadlocks 25 It is assumed that the technological branches have the same probability, so OGAnalyzer assigns the reciprocal of the number of branches to these edges. The fault caused branches are collected into a table and the user has to define the probability of faulty and normal modes as it can be seen in Fig.11. The occurrence graph with arc weights depicted in OGAnalyzer can be seen in Fig.12. The black numbers attached to the arcs are the arc weights. Values equal to 1 were not depicted. Using these probability values the OGAnalyzer calculates automatically the probability of a node by multiplying arc weights along the path leading from the node representing the initial state to that node. If two or more paths lead to the node then the probability values belonging to these paths are summed. The calculated probability values can be also seen in Fig.12 as red numbers assigned to nodes. In this manner the probability of all system states in the investigated system can be calculated. For example, the probability of faultless completing of both workpieces is 0.28 + 0.28 = 0.56, which is equal to the sum of the probabilities belonging to the nodes 55 and 57. The probability of that case when the first workpiece is manufactured without fault, but the label of the other piece gets damaged during the second manufacturing process is 0.03 + 0.03 = 0.06 (sum of the probabilities belonging to the nodes 46 and 48). In case of timed nets the occurrence graph is often a tree and the same token distribution belongs to different nodes because of the different time stamp. For the determination of probability of a given system state, the probability values belonging to different nodes have to be summed up. Figure 11: Defining the probabilities values of faults Figure 12: The occurrence graph with arc weights and probability values Figure 10: The occurrence graph of the example in the window of OGAnalyzer 26 Conclusion A novel occurrence graph investigation procedure for discrete event systems described by Petri nets was proposed in this paper for model-based diagnostic purposes that utilize the knowledge of the occurrence probability of faults. The model of the investigated system was defined in timed coloured Petri net form. The colours of tokens representing the workpieces were used to distinguish them and to assign a label of the processes to be carried out. The arc inscriptions and built-in probability functions were used for the fault modelling and integration. The operation of the system was investigated via simulation both in timed and non-timed cases and in faultless and possible fault operational modes with given number pieces and prescribed manufacturing lines. For the behavioural analysis of the model the occurrence graph method was used. A special software module, called OGAnalyzer has been developed for the handling of the probabilities on the occurrence graph and for calculating the occurrence probability of system states. Acknowledgements This research is partially supported by the Hungarian Research Fund through grant No. K-83440. We also acknowledge the financial support of the Hungarian State and the European Union under the TAMOP- 4.2.2.A- 11/1/ KONV-2012-0072. REFERENCES [1] CASSANDRAS C.G., LAFORTUNE S.: Introduction to Discrete Event Systems, Kluwer Academic Publishers, 1999 [2] GERZSON M., MÁRCZI B., LEITOLD A.: Diagnosis of Technological Systems based on their Coloured Petri Net Model, ARGESIM Report no. S38 (Eds. TROCH I., BREITENECKER F.) 2012, p. 358/1–6 [3] JENSEN K.: Coloured Petri Nets: Basic Concepts, Analysis Methods and Practical Use, Springer- Verlag, 1997 [4] CPN GROUP, University of Aarhus, Denmark: CPNTools 2.2.0 http://wiki.daimi.au.dk/cpntools/ (last accessed: May 25, 2014) [5] FANTI M.P., SEATZU C.: Fault diagnosis and identification of discrete event systems using Petri nets, Proc. 9th International Workshop on Discrete Event Systems, WODES, 2008, 432–435 [6] JENSEN K., KRISTENSEN L.M., WELLS L.: Coloured Petri Nets and CPN Tools for Modelling and Validation of Concurrent Systems, Int. J. of Software Tools for Technology Transfer, 2007, 9(3–4), 213–254 [7] MURATA T.: Petri Nets: Properties, Analysis and Applications, Proceedings of the IEEE, 1989, 77(4), 541–580 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 27–34 (2013) IMMEDIATE EVENT-AWARE MODEL AND ALGORITHM OF A GENERAL SCHEDULER TIBOR DULAI1!, ÁGNES WERNER-STARK1, AND KATALIN M. HANGOS1,2 1 Department of Electrical Engineering and Information Systems, University of Pannonia, Egyetem str. 10, Veszprém, H-8200, HUNGARY 2 Process Control Research Group, Computer and Automation Research Institute, Kende u. 13-17, Budapest, H-1111 HUNGARY !E-mail: dulai.tibor@virt.uni-pannon.hu A stochastic scheduling problem is investigated in this work that considers workpieces to be manufactured according to individual recipes containing manufacturing steps performed by workstations as resources. Unexpected stochastic breakdown of a workstation or the faulty termination of a recipe, when a manufacturing failure renders the workpiece out of specifications, forms the set of immediate events. A model and an algorithm are proposed as the basis of a scheduler, which takes into account the possible immediate events, estimates their probability and suggests resource allocations which provide the best overall work-flow even when an immediate event happens. This model includes the possibility of handling alternative resources that can substitute each other in case of an immediate event, like sudden technical failure of a resource. Immediate events are not exactly predictable; however, based on previous experiences, their probabilities can be estimated. Our model uses the properties of the resources (including how they can substitute other types of resources) and the required sequence of them during the workflow (i.e. the recipes). The proposed scheduling algorithm constructs a solution workflow that reacts in the best way (in average) even for an unexpected event. The proposed model and scheduling algorithm is illustrated on two industrial case studies. Keywords: scheduling, resource allocation, alternative resources, immediate event-awareness Introduction Scheduling is an important and widely used topic of operations research. Besides of its theoretical importance, industries can also benefit from optimal schedules and resource allocation. Several different algorithms were established for organizing process elements on the time scale [1, 2], related to for example computer networks [3], business processes [4, 5] or industrial processes [6]. The results are usually represented on a Gantt chart. There are also differences between these algorithms related to their application times. Some of them are applied offline before the scheduled processes starts, others are applied real-time. Real time scheduling techniques, and a special adaptive real-time scheduler is introduced in Ref. [10]. For improving their results, there are cases when the methods use historical data during the creation of the schedule [11]. The common scheduling methods usually handle the resources individually and do not take into account the relationships, e.g. the similar functionality between them. Only few publications investigate the cooperation possibilities that are enabled by the similar functionality of resources in scheduling [7]. In this paper a scheduling method is proposed, where resources may substitute each other and immediate failures of resources may happen. After we introduce the problem and its main building blocks, we present our model and the algorithm developed. Its operation and properties are demonstrated on a simple and a more complex problem as case studies. Problem specification A general scheduler intends to determine the placement of activities of resources on the time scale. In an advanced case one might consider additional features given in the problem’s model – like substitution possibility of the resources – which help to redefine the classic scheduling problem to be usable in different real- life applications. This substitution may help in achieving a certain fault tolerance property in such a way that a technical failure causes the least possible negative effect. The basic aspects and sub-tasks related to this extended scheduling problem are collected in this section. 28 Cooperation In some cases more than one resource is able to carry out a particular activity (usually with different productivity). It means that when the appropriate resource is busy, there can be another resource that is able to take over a particular task if required. We shall term this case a cooperative situation when substitution is possible, In other words, the resources can cooperate with each other. Taking the cooperation aspect into account during scheduling the performance may increase [8]. We introduced the so called substitution vector as an element of our model to handle substitution related sub-problems. Any scheduler [1, 2] can be applied as a basis of the cooperative extension, which inputs the sequence of resources as input processes and has to arrange them on the time scale taking into account the availability and temporal constraints. All the improvements are done on this basic schedule. Improving fault tolerance When a schedule is ready and the workflow starts, an immediate event (e.g. a technical failure of a resource) in a process may influence other processes, too, in a negative way [9]. If the effects of that event are calculated in the schedule, then the solution with the best answer to immediate events can be selected from the set of solutions with the same performance, and the faults’ negative impact on the schedule decreases. Our model deals also with this aspect. The model For the above mentioned problem set, we created an universal model and an algorithm to carry out the scheduling. The work intends to be the basis of several different tasks with the goal that the operation time should be minimal. In this section we introduce our model and the scheduler that works on this model. Main parameters, notations and functions The model and its parameters are designed in such a way that it can be applied for different kinds of scheduling problems (e.g. scheduling of industrial production processes, test processes, scheduling and managing the resources of electrical networks, etc.). For this reason we collected the necessary parameters that make possible the development of a general framework, which uses different kinds of resources in different processes. It is able to handle different needs. In this sub-section we introduce the main parameters and functions we use in our work. These are as follows: Process = {proc1, …, procj}, (1) stands for the set of processes; P = {p1, …, pl}, (2) stands for the set of product types; R = {r1, …, rm}, (3) stands for the set of resources; O(ri) = {oi1, …, oik}, (4) stands for the set of the operation modes of resource i; A = {a1, …, an, pause(t)}, (5) stands for the set of basic activities, where pause(t) is an empty activity with length of t hour; ra : R × A → {0, 1}, (6) defines a function for determining whether a resource is able to perform an activity; rap: R × A × P → {0, 1}, (7) defines a function for determining whether a special activity of a resource can be applied for a product type; t: R × A × P × O → N, (8) provides the suggested operation time of a resource in a given operation mode performing a given activity on a given product type; q: R × A × P × O → [0,100], (9) provides capacity information: how many percent of a resource capacity is occupied by one piece of a given product type in a given operation mode of a resource while performing a given activity; e: R × A × P × O → [0,100], (10) provides the probability of resource failure during performing a given activity on a given product type in a given operation mode; s: R × A × P × O → {0, 1}, (11) results in a binary decision: whether a resource activity in a given operation mode on a given product type can be suspended without restarting it from its beginning; rreq: R × A × P × O × DateTime → P(R, N) (12) (power set on pairs of a resource and a natural number), provides the additional resource need of a resource’s given activity on a given product type in a given operation mode in a given hour subst(ri,aj, pk,oil ) = [ hi1 ... him ] , (13) 29 is the substitution vector, where m is the number of resources, hin is a natural numbers for all 1 ≤ n ≤ m, ri is the ith resource, aj is a basic activity, pk is a product type, oil is an operation mode of resource i, hin denotes how many percent of the productivity of resource i is needed in a given activity in a given operation mode on a given product type to substitute totally resource n supposing unchained operation time. The substitution vector is calculated based on the function ra and other functions (e.g. function t in case of time-based optimization). map: R x DateTime → N, (14) provides the information on accessibility of resources: how many resource of a given type is accessible in a given hour; F: R → N, (15) provides the expected number of hours how long a given resource is unavailable in case of its failure; tstart: Process → DateTime, (16) shows the start time of a process; tmaxend: Process → DateTime, (17) shows the maximum finish time of a process; prev: Process × A × N+ → {A, Ø}, (18) provides the prior basic activity of the nth occurrence of a given basic activity in a given process; next: Process × A × N+ → {A, Ø}, (19) provides the following basic activity of the nth occurrence of a given basic activity in a given process; maxdelay: Process × A × N+ → N, (20) provides the maximal duration of time out, which is tolerated by the nth occurrence of a given basic activity in a given process; dur: Process × A × N+ → Q, (21) provides the time scale, which number the default operation time of the of the nth occurrence of a given basic activity in a given process has to be multiplied with, for getting the real operation time of the activity. The proposed algorithm Our algorithm was created to take into account the cooperation possibilities of the resources and to have failure-aware behaviour during scheduling. The algorithm can be separated into three main blocks, as it is shown in Fig.1. In the first phase, we associate resources for each basic activity of the processes. At the start, processes are described only as a sequence of basic activities. We have to turn it into the “language” of resources. In this task, we use the function ra for determining, which resources are able to carry out the desired basic activity. Furthermore, we deal with the minimization of the operation time; we select the resource, which has the operation mode with the minimal operation time, i.e. ri with min(t(ri, A, P, oil )), ∀oil ∈O(ri ), ra(ri,A) =1 (22) is selected. The first step of our algorithm results in a sequence of resources for each process. The second phase of the algorithm is performed by the main scheduler task, taking into account the alternative resources. The first task in this phase is to select the “basic” process, which has the least robustness. We do it by selecting proci with min(tmaxend (proci )− tstart (proci )−TDUR) , (23) where TDUR is the sum of the durations of all of the basic activities of proci. After the determination of the basic process, the algorithm enters a loop, in which it selects a process from the set of the remaining processes and attempts to place all the basic activities of the selected process element-by-element on the time scale. During this operation the algorithm handles the substitutability of the resources. The insertion of a basic activity starts at the initial activity of the process and each activity is inserted into the earliest possible time point. This assumption is a fundamental point in the method. The insertion of an activity at its earliest time may have two possible outcomes: • successful insertion • unsuccessful insertion, which means that there is not enough amount of the resource to serve the activity which starts at its earliest time point. We call this case a collision. In case of a collision, the algorithm attempts to find an alternative resource using the substitution vector, Figure 1: The basic parts of our algorithm 30 which has enough leftover capacity to carry out the activity at the desired earliest time point. If this is not possible, because there is not enough resource at the desired time, the algorithm tries to shift one of the collided activities to solve the problem. The detailed mechanism of the activity-insertion task is illustrated in Fig.2. When the activity-shift resulted in no success, the algorithm attempts to repeat the same process with other alternative resources. If there is no success, the whole process will be left from the final solution; otherwise the initial scheduling is ready. In the final, third phase of our algorithm the schedule is tried to be modified in such a way that it should have a higher fault tolerance than the original one. This phase assumes that fault happens at the processing of special activities of the processes and calculates its effect into the schedule. We do not want to hurt the time constraints of the original problem that’s why one of the main parameters of this phase are the maximum finish time values (tmaxend) of the processes. We process the processes starting from their end in the following way: Let WP be the process list we work on. At the start WP = Process. For each element of WP, we create a pointer, which points at the penultimate activity of the process. The activities pointed by the pointers (let their set be WA) are candidates for assuming them to have fault. Starting from this initial state, this phase of the algorithm works as follows: 1. Select the aj activity from WA with the highest error probability: max(e(ri, a j, pk, oil )), a j ∈WA (24) 2. If the start of the basic activity which follows ai can be shifted by ai’s duration plus by F(R(ai)), where R(ai) is the resource which carries out ai without causing collision and without exceeding tmaxend(prock), where prock is the process which belongs ai to, then (i) the shift will be done, and (ii) the basic activity, which precedes ai in the process will be added to WA and ai will be taken out from there. It means that the pointer of the process will be set one step backward. Else, ai is taken out from WA and its process is taken out from WP, too. It means that the activities of the process are not modified any more. This phase of the algorithm modifies the placement of the basic activities of the processes on the time scale in such a way that there is no need to change the scheduled activities in case of failures with high probability. It means that these failures will not cause delay and exceed of the maximum finish time. Case studies In this section, the previously presented algorithm is applied to two hypothetical manufacturing examples: we show how it takes into account the substitutability of the resources and how to make the schedule ready for possible technical failures. A simple problem and its solution Let there be four processes with the same start and maximum finish time (tstart = 0 and tmaxend = 65 are identical in each cases). The sequences of the process activities can be seen in Fig.3. In this example, we work only with one product type (P), and we intend to minimize the maximum operation time. Let’s suppose that in this example there are two pieces of resource r1, while only 1 piece of the others, and ra(ri, aj) = 1 only if i=j, (25) excepting two cases: ra(r2, a3) = 1 (26) and ra(r3, a2) = 1, (27) as well. Moreover, the substitution vector of r2 is: subst(r2,a3,P,O) = [ 0 100 125 0 ] , (28a) while the substitution vector of r3 is: subst(r3,a2,P,O) = [ 0 200 100 0 ] , (28b) Figure 3: Activity sequence of the example processes Figure 2: Activity insertion onto the time scale 31 which means that r2 can substitute r3, and r3 can substitute r2. However, the substitution results in 125 percent or 200 percent more time, respectively. Another assumption is when only r1 has two operation modes: in one of its operation mode its operation time at a1 activity on product P can be reduced with 10%. As our goal is to minimize the necessary time, we will use this operation mode of r1 in each case. Moreover, assume that all the resources are unavailable for 2 time units in case of failure. Based on this information, the first step of our algorithm creates the sequence of the resources for each process. This sequence completed with the durations of each activity after user-based modification of the default operation times is illustrated in Fig. 4. The second phase of the algorithm determines the basic process; in our case it is proc3 with its 63 unit length. Its representation on the time scale is shown in Fig.5. After this selection a loop is started, choosing the less robust process from the remaining set. The next process is proc1. Its first activity can be placed onto the time scale without collision. The placement of the second activity collides with the basic process’s first activity. The first reaction of the algorithm is to search for alternative resource. Resource r2 can substitute r3 without collision; however in this case the operation time will increase from 4 to 5 units. The third activity of proc1 can be inserted, since there are two pieces of resource r1. The fourth activity collides again; however substitution solves this problem, and the problem-free placement of the final activity of the process can be seen in Fig.6. The third process to handle is proc4. Its first activity can be placed easily; however, the second activity has collision and there is no possibility to substitute it without collision either. This case requires the shift of the activity as illustrated by Fig.2. The activity that shifts the later start is reason for a pause was defined in proc4. The third and fourth activities have to be also shifted. The difference between the two cases is the resource that is able to substitute gets free earlier. This is why the fourth activity will be carried out by r2 instead of r1, started after a short pause. After the placement of the process’s final activity and all the activities of proc2,, this provided the schedule with some activities of the previously placed processes shifted necessarily twice, as illustrated in Fig.7. The final phase of the algorithm is to make the schedule to be fault-aware. As tmaxend = 65 and the length of proc3 is 64, we will not modify this process. Similarly, we will not modify proc1 and proc4 either. Only proc2 lets the algorithm to prepare it to be fault- aware, and the modification can be applied to all of its Figure 5: Placement of the basic process onto the time scale Figure 6: The schedule after handling proc1 Figure 7: The schedule after the second phase of the algorithm Figure 8: The final schedule Figure 4: Sequence of the initial resource allocation 32 activities. The result can be seen in Fig.8. A more complex problem and its solution In this example we deal with six processes. In contrast with the previous example, they differ in their start and maximum finish time as shown in Table 1, where the values are represented in time units. The activity- sequences of the processes are illustrated in Fig.9. In the example form Fig.9, eight different activities are applied. The next question is which resources are able to carry out these activities. In this example, there are six kinds of resources, six resource types (r1 - r6). Suppose that two pieces of resource type 2, 3 and 5 exist, while all the other resource types have only one representative. The substitution vectors of the resources are as follows: subst(r1,a2,P,O) = [ 100 150 0 0 0 0 ] , (29a) subst(r2,a1,P,O) = 150 100 0 0 0 0! " # $ , (29b) subst(r3,a5,P,O) = 0 0 100 0 300 0! " # $ , (29c) subst(r5,a3,P,O) = 0 0 120 0 100 0! " # $ , (29d) subst(r4,a1,P,O) = 166 0 0 100 0 0! " # $ . (29e) It means that r1 and r2 resources may substitute each other, r3 and r5 resources are able to substitute each other and r4 resource is capable to substitute r1 resource. Other substitutions are not possible. These vectors show that ra(r2, a1) = 1, (30) ra(r1, a2) = 1, (31) ra(r3, a5) = 1, (32) ra(r5, a3) = 1 and (33) ra(r4, a1) = 1. (34) Moreover, we suppose that ra(ri, aj) = 1, if i = j; (35) ra(r2, a7) = 1 and (36) ra(r5, a8) = 1. (37) This means that a7 activity can be carried out by r2 resource, and a8 activity can be carried out only by r5 resource. In case of failure, the resources are unavailable as much as shown in Table 2. We suppose that each resource works only in one operation mode and all of the resources have the same error probability. In the example we work only with one product type (P), and we intend to minimize the maximum operation time. Applying the first step of our algorithm, we create the sequence of the resources for each process. This sequence completed with the durations of each activity after user-based modification of the default operation times as illustrated in Fig.10. At this phase we apply the second step of our algorithm and determine the robustness of each processes. In the calculation, we use the following computation method, as we mentioned earlier: tmaxend (proci )− tstart (proci )−TDUR (38) The values obtained are summarized in Table 3. Table 1: Some properties of the processes of the second case study process start time, tstart max. finish time, tmaxend p1 0 30 p2 15 30 p3 7 30 p4 3 25 p5 0 23 p6 6 29 Figure 9: The activities of the second example’s processes Figure 10: Resource-sequence of the example’s processes Table 2: The effect of failure on the second example’s resources Resource r1 r2 r3 r4 r5 r6 Time of unavailability (F) [time unit] 2 1 3 2 2 1 33 Table 3 indicates that the basic process is proc4, which is the least robust process. We place its activities onto the time scale as a chart illustrates this in Fig.11. The black line represents the 0 time point and the red one signs the maximal finish time of the process. The second least robust process is proc3. The placement of its activities onto the time scale can be done easily, without any collision as shown in Fig.12. The third least robust process is proc2. During its placement there is one collision: its second activity collides with proc4’s fifth activity. However, it doesn’t cause any problem, because there are two pieces of resource r2. The second resource of r2 resource type can carry out the activity in the originally planned time. The result is shown in Fig.13. Looking for the next least robust process, proc6 is the next. Its fourth activity collides, thus the second piece of r3 resource type has to be used for carrying it out. Moreover, its final activity also collides. Unfortunately, there is only one piece of r1 resource type. However, r2 resource type can substitute an r1 resource. The problem is that both r2 resources are occupied at the desired time. There is another resource type, r4, which is able to substitute an r1 resource with a bit more necessary time than r2. As r4 is free at the desired time, it will take over the task. The time it requires for its task is 1.66 times more than an original r1-type resource would need for that. Fig.14 illustrates the results. The second most robust process is proc5. During the placement of its activities onto the time scale, we find that the second activity collides with proc4’s second activity. Fortunately, there is a second piece of r3 resource type which can be used freely for the desired time interval. It solves the problem and results a chart as seen in Fig.15. At last, three problematic cases happen during process 1: the collision of the third and fourth activities can be solved by the second piece of the desired resource types; however, in case of the final activity only time shift can solve the collision, because there is only one piece of r4 resource type and it can not be substituted by any other resource. The result is shown in Fig.16. After placing all activities onto the time scale and solving all collisions, only the fault-tolerance-related improvement need to be done, as the final step of our algorithm. This phase intends to shift the activities of the processes in time, starting from their end by the duration of their preceding activity plus its F value. Because of the maximal finish time constraint, proc4 and proc3 (the first two processes on Fig.16) cannot be modified. If we look at Fig.16’s third process, its last activity requires F(r1) = 2 time units plus its normal operation time in case of failure. If we take into account the maximal finish time of the process, its place cannot be modified. This is the case for the fourth process, as well. The fifth process in in Fig.16 ends with r3 resource. In case of a failure, it requires 8 time units until it reaches the maximal finish time of the process. Taking this into account, only 3 time units remain, Figure 11: Placement of the second example’s basic process onto the time scale (proc4) Figure 12: Placement of the second example’s second process onto the time scale (proc3) Figure 13: Placement of the 2nd example’s third process onto the time scale (proc2) Figure 14: Placement of the second example’s fourth process Figure 15: Placement of the second example’s fifth process onto the time scale (proc5) Figure 16: Placement of all processes of the second example onto the time scale (proc1) Table 3: The robustness values of the example’s processes process proc1 proc2 proc3 proc4 proc5 proc6 robustness 16 4 2 0 11 8 34 which is not enough for the previous activity. That’s why we do not modify this process, either. The final activity of the last process, proc1 is carried out by an r4 resource. In case of a failure, it requires 4 extra time units. In order to avoid collision we shift this activity a bit more forward as shown in Fig.17. As F(r2) = 1, and a2 activity requires 3 time units to be carried out by r2 resource, proc1’s fourth activity can also be shifted in time. Its previous activity would require 6 time units after itself. Its shift would result 1 extra time unit for previous activity, which is not enough for making it ready for tolerating fault. This hinders building in more fault tolerance. The final result of the algorithm is shown in Fig.18, in which we have a schedule, which took into account the possible substitutions, the capabilities of the resources; moreover in some places, it tolerates faults of resources without the need of rescheduling. Conclusion We presented a model and algorithm for creating schedules, which tolerate some resource-failures and utilize the substitutability possibilities of the resources. Our goal was to establish a model, which can be the basis of applications in different segments of scheduling problems and makes it possible to generate schedules with the optimization criteria related to the operation time. We introduced an algorithm, which creates an initial schedule on the basis of the input processes and the substitutability of the resources. This can be improved with preparation for likely failures of resources. The operation of the algorithm was presented on two hypothetical case studies. Acknowledgements This research was supported by the European Union and Hungary and co-financed by the European Social Fund through the project TÁMOP-4.2.2.C-11/1/KONV-2012- 0004 and TAMOP-4.2.2.A-11/1/KONV-2012-0072, and in part by the Hungarian Research Fund (OTKA) grant number K83440. Support from the National Research Centre for Development and Market Introduction of Advanced Information and Communication Technologies is also acknowledged. REFERENCES [1] PINEDO M.L.: Scheduling: Theory, Algorithms, and Systems, 4th ed., Springer, New York, USA, 2012 [2] BRUCKER P.: Scheduling algorithms, 5th ed., Springer, Osnabrück, Germany, 2007 [3] BRAHIMI B., AUBRUN C., RONDEAU E.: Modelling and Simulation of Scheduling Policies Implemented in Ethernet Switch by Using Coloured Petri Nets, Proc. 11th IEEE Int. Conf. Emerging Technologies and Factory Automation, Prague, Czech Republic, 2006, 667–674 [4] BARBA I., VALLE C. DEL: A Job-Shop Scheduling Model of Software Development Planning for Constraint-based Local Search, Int. Journal of Software Engineering and Its Applications, 2010, 4(4), 1–16 [5] XU J., LIU C., ZHAO X., YONGCAREON S.: Business process scheduling with resource availability constraints, Proc. OTM’10, Hersonissos, Crete, Greece, 2010, 1, 419–427 [6] SULE D.R.: Production Planning and Industrial Scheduling: Examples, Case Studies and Applications, 2nd Ed., CRC Press, Ruston, LA, USA, 2007 [7] MURTHY S., AKKIRAJU R., RACHLIN J., WU F.: Agent-based cooperative scheduling, Proc. AAAI Workshop on Constraints and Agents, Providence, RI, USA, 1997, 112–117 [8] DULAI T., WERNER-STARK Á.: Immediate event- aware routing based on cooperative agents, Proc. Factory Automation, Veszprém, Hungary, 2012, 144–148 [9] PALOMBARINI J., MARTINEZ E.: SmartGantt – An intelligent system for real time rescheduling based on relational reinforcement learning, Expert systems with Applications: An International Journal, 2012, 38(11), 10251–10268 [10] NANDANWAR J., SHRAWANKAR U.: An adaptive real time task scheduler, IJCSI International Journal of Computer Science Issues, 2012, 9(6/1), 335–340 [11] GREGG C., BOYER M., HAZELWOOD K., SKADRON K.: Dynamic heterogeneous scheduling decisions using historical runtime data, Proc. 2nd Workshop on Applications for Multi- and Many-Core Processors, San Jose, CA, 2011 Figure 17: Building fault tolerance into the last process of the second example Figure 18: Building fault tolerance into the last process of the second example HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 35–40 (2013) ESTIMATION OF PARAMETERS FOR AN EXTREMELY LOW FUEL CONSUMPTION INTERNAL COMBUSTION ENGINE-BASED MEGAMETER-III VEHICLE ISTVÁN PINTÉR1! AND MIHÁLY BAGÁNY2 1 Department of Informatics, Kecskemét College, Izsáki str 10, Kecskemét, H-6000, HUNGARY 2 Department of Natural Sciences and Engineering, Kecskemét College, Izsáki str 10, Kecskemét, H-6000, HUNGARY !E-mail: pinter.istvan@gamf.kefo.hu Megameter-III is a special vehicle designed and constructed for international competitions by our team. A crucial part of the vehicle is the internal combustion engine, which ought to have extremely low fuel consumption. The experiences of previous competitions confirmed that it is necessary to develop a professional telemetry system, which allows tracking and logging online important parameters of the engine and the vehicle. In this paper, we discuss results for estimating of several parameters of the engine (angular velocity, angular acceleration, the Coulomb-Morin external friction torque and the damping constant). In our work, injector-pulses and inductive sensor’s signals were monitored. Keywords: internal combustion engine, extremely low fuel consumption, angular velocity, angular acceleration, external friction torque, damping constant Introduction Since 2010, the racing team at Kecskemét College participates in international competitions of fuel- efficient vehicles. The vehicles were named after 1 Megameter being 1000 km. The first result with Megameter-I in gasoline category was 1588 km L-1 (2010, Lausitz, Germany), and the best result of Megameter-II was 2661 km L-1 (2011, Nokia, Finland). Megameter-III achieved 2696 km L-1 (2012, Rotterdam, Nederland. The team-photo can be seen in Fig.1). The latest model, the Megameter-IV has exceeded the 3000 km/litre dream-limit (3082 km L-1; 2013, Nokia, Finland) [1]. During the development of the Megameter series, it became clear that in order to improve the technical parameters of the vehicle it is necessary to establish a telemetry system suitable for measuring, displaying and logging the engine’s and vehicle’s parameters during the race. This paper presents our recent results in this work, which have already been utilized in developing and construction of the Megameter-IV. The measurements were taken using the Megameter- III engine, which has the following characteristics: • one cylinder, four-stroke, air-cooled • overhead controlled, overhead valve, intake manifold injection • bore/stroke: 31.5/45 mm • displacement volume: 45 cm3 • compression ratio: 14 • power: 700 W • torque: 2.2 N·m • specific fuel consumption 235 g kWh-1 There are different application areas and methods for estimation of engines’ angular velocity. These estimation procedures are based mainly on crankshaft’s rotation angle [2, 3, 4]. In the problem of model-based estimation of mechanical losses, the dependence of instantaneous angular velocity on angular position plays an important role [4]. A second-order spline- interpolation-based method was proposed in [2] for real- time measurement of angular velocity and angular acceleration. Our telemetry system is related to this latter work, and our results proved to be useful in further developing our special vehicle. Figure 1: Student team photo with the Megameter-III vehicle (Rotterdam, 2010, second place) 36 The physical model with its approximations and the parameter estimation In order to achieve ever decreasing fuel consumption it is necessary to uncover the sources of energy losses. For this purpose a suitable physical model of internal combustion engine (ICE) is needed [5]. Moreover, with an appropriate physical model in hand, one can design the measurement for estimating the parameters of the physical model. Physical model for estimating the external friction torque and damping constant The instantaneous angular velocity of the crankshaft is changing even during only one crankshaft revolution. It is smaller on the compression stroke and greater in the midst of power stroke. The friction torque affecting the engine, averaged for two crankshaft revolutions depends on the instantaneous angular velocity. By expanding this function up to the second order we get Eq.(1): −M ≈ M 0 + M1 ω1 ω + M 2 ω2 2 ω 2 (1) where the indexed quantities are constant and positive, and they have physical meanings as follows: M0 is the Coulomb-Morin external friction torque, which is independent of the relative speed of frictional surfaces (in this case the relative speed is the angular velocity); M1/ω1⋅ω defines the damping constant (the second term corresponds to torque resulting from Newton’s internal friction’s force, it is proportional to the first power of angular velocity); M2/ω2 2⋅ω2 stands for the drag torque characterizing the turbulent flows in pipes; it is in proportion to the square of angular velocity and it is resulting from the flow resistance in the intake and exhausting channels, in the valve-cage and in the air- pipe of the crankcase. Our task is the estimation of the external torque from the measured data; therefore it is necessary to estimate the angular velocity. The basic equation of rotational motion is shown in Eq.(2): −M = −Θ dω dt (2) where after substitution of Eq.(2) into Eq.(1) we get M 0 + M1 ω1 ω + M 2 ω2 2 ω 2 = −Θ dω dt (3) Let’s assume that the initial value of the angular velocity is ω0. By separating the variables and integrating we get Eq.(4): dω M 0 + M1 ω1 ω + M 2 ω2 2 ω 2ω0 ω ∫ = − dt Θ0 t ∫ . (4) The result of the integration is shown in Eq.(5). t = − 2Θ 4ac− b2 atan 2cω + b 4ac− b2 # $ % & ' (+ + 2Θ 4ac− b2 atan 2cω0 + b 4ac− b2 # $ % & ' ( , (5) where three new variables were introduced for more readability as shown in Eq.(6). a =M 0 , b = M1 ω1 , c = M 2 ω2 2 2 . (6) By solving Eq.(5) for ω one can get the final solution. However, instead of doing so, we analysed two approximate solutions of Eq.(4). As we will see below, these approximate solutions fit real measurements well. In order to estimate the engine’s friction torque and damping constant long measurement duration (e.g. 20 s) is necessary. The relevant physical quantities can be estimated from the measured injector pulses. Approximation 1: the friction torque is independent from angular velocity By assuming that the other two terms are negligible compared to external friction M0 the differential Eq.(4) simplifies to: dω M 0ω0 ω ∫ ≈ − dt Θ0 t ∫ , (7) and its solution is: ω =ω0 − M 0 Θ t , (8) where ω0 is the value of angular velocity at the beginning of deceleration. The angular velocity of the engine’s crankshaft has to be measured as the function of time for parameter estimation. Namely, after line-fitting to measured data the intercept and slope can be computed from relationship ω̂ = A− Bt , so both the estimation of Coulomb-Morin external friction torque and the estimation of initial value of angular velocity can be determined: ω0 = A , M 0 = BΘ . 37 Approximation 2: the friction torque is independent from the squared value of angular velocity When the engine is running with wide open throttle, and the friction torque of pipe-flow losses is negligible comparing to external and internal losses, the differential equation to be solved is the following: dω M 0 + M1 ω1 ωω0 ω ∫ ≈ − dt Θ0 t ∫ . (9) The solution of Eq.(9) is: ω = ω1 M 0 M1 +ω0 ! " # $ % &e − t τ −ω1 M 0 M1 , (10) where τ = ω1Θ/M1 is the time constant characterising the angular deceleration. By measuring the engine’s crankshaft’s angular velocity as the function of time, and after fitting a biased exponential-type curve we get the estimation CAe Bt −= −ω̂ , from which the estimates of physical quantities are the following: - the external friction torque is: M0 = BCΘ, - the damping constant is: M1/ω1 = BΘ, - the initial value of angular velocity resulted from curve-fitting is: ω0 = A-C. Estimation of the instantaneous torque and the running irregularity The instantaneous angular acceleration as the derivative of the instantaneous angular velocity is needed for estimation of the instantaneous torque. Moreover, knowing the instantaneous angular velocity, another important parameter of the engine, the running irregularity can be determined. From the point of view of measurement, the estimation of instantaneous angular velocity requires the duration of two crankshaft revolutions (e.g. 100 ms). For measuring this, an inductive sensor was used. When determining the running irregularity, the engine and the dynamometer were disconnected and the flywheel ring gear teeth together with an inductive sensor served as signal source. The induced voltage from the sensor can be approximated as: U U0 = ω +β t( ) ω0 cos zω t + z β 2 t 2 ! " # $ % & , (11) where U is the output signal of the sensor, U0 is the amplitude of the sensor corresponding to the constant angular velocity of ω0, ω is the angular velocity of the crankshaft, β is the angular acceleration of the latter, z is the number of teeth of ring gear. In case of zero or small values of angular velocity the relationship Eq.(11) simplifies to U U0 = ω ω0 cos zω t( ) . (12) By assuming that during the pass-time of one tooth (tooth-time) the angular velocity is not changing significantly, the signal amplitude during the tooth-time can be well approximated by a constant. That is, the approximated sensor signal when passing the zth tooth is given as U/U0k = cos (z ωk t). The running irregularity is defined as a ratio in Eq.(13): δ = ωmax −ωmin ωaverage (13) Measurements, estimations and evaluations The measurements were accomplished in the Student’s Workshop of GAMF Faculty using an Agilent DSO-X 2002A digital storage oscilloscope. The oscilloscope is two-channel type with upper input cut-off frequency of 70 MHz, maximum sampling frequency of 2 GS/s and with sample storage of 50000 samples. The algorithms used for evaluation of measured data were intentionally different. In one case the data series were not re-sampled and the differentiation was approximated by finite differences, and in the other case both re-sampling and smoothing derivative algorithm were applied. Moreover, the computation method of the trend of angular acceleration was also different. Estimation of the external friction torque and damping constant In this setup, an extra flywheel was driven by the engine fixed on the test stand. The measurement setup can be seen on Fig.2. The reduced resultant moment of inertia of synchronously rotating components with the crankshaft was Θ ≈ 0.0534 kg m2. The dynamometer’s friction torque was approximately constant with a value of Mdynamometer ≈ 0.2 N m. Figure 2: Arrangement for measuring the friction torque 38 The engine’s injector was triggered by a signal from the crankshaft, one injection occurred in every two crankshaft revolutions. The oscilloscope was connected to the injector’s solenoid, and the solenoid’s signal was sampled with a sampling frequency of 2.5 kHz, that is the sampling interval was Ts = 400 µs. The part of measured injector pulse-signal can be seen on Fig.3. The crankshaft’s angular velocity ω = ω (t) as a function of time was determined from the stored samples, with duration of 20 seconds. At the beginning of the measurement, the engine was accelerated to 4000–4200 rpm, then the injector was disconnected from the intake manifold, and the engine was allowed to decelerate freely. By defining a suitable trigger level, the elapsed time T between two consecutive rising edges of the impulse can be determined. The instantaneous frequency is two times the reciprocal value of T (finstantaneous = 2 / T), because every two revolution results in one impulse. The instantaneous angular velocity can be computed from this instantaneous frequency value. As the resulting angular velocity value sequence is non- equidistant, it has been re-sampled using cubic spline interpolation. First, a linear model was fitted (Fig.4). The results are summarised in Table 1. The two methods gave similar parameters, but the resampling resulted in somewhat better goodness-of-fit. The estimation of the physical quantities from parameters shown in Table 1 can be seen in Table 2. The estimations resulted by the two methods are similar in case of model of Eq.(7). In case of model of Eq.(9), the fitting has been fulfilled on the resampled data series. The results are demonstrated on Fig.5. From the above results the external friction torque is M0 = B⋅C⋅Θ ≈ 0.419 N·m; the average braking torque of the engine is Mengine = M0 - Mdynamometer ≈ 0.219 N·m; the damping constant is M1/ω1 = BΘ = 6.94⋅10-4 N·m·s. The initial angular velocity and revolution are ω0 = A – C ≈ 352 rad s-1, n0 = 3358 rpm. We can conclude that the two types of models gave similar results, but the biased exponential model fitted somewhat better the measured data than the linear model, thus the former are preferred. Furthermore, the presented data indicate that both the average breaking torque of the engine and the damping constant are acceptable small, therefore the engine is acceptable from this point of view. Estimation of the instantaneous torque and the running irregularity from measurements When the running irregularity was measured, the engine and the dynamometer were disconnected. The moment of inertia of the parts running synchronously with the crankshaft was Θ = 4000 kg·mm2. A ring gear fixed to the crankshaft was used as signal source (number of teeth is z = 66) together with an inductive sensor fitted from the teeth in a distance of 4 mm. An oscilloscope was connected to the sensor and the voltage signal was sampled with 500 kHz sampling frequency (Fig.6). It is necessary to use such a great sampling rate because of the following. In order to acceptably estimate the instantaneous angular velocity, it is necessary to measure N samples during the tooth-time. In case of nominal revolution of n0, the sampling frequency is fS = Figure 3: Measured injector pulses using a power MOSFET Table 1: Results of the fitting of data shown in Fig.4 Resampling A, rad s-1 B, rad s-2 R2 No 353.78 12.103 0.9954 Yes 346.75 10.952 0.9981 Table 2: Estimation of physical quantities Resampling n0, rpm M0, N·m MICE, N·m No 3378 0.650 0.450 Yes 3311 0.585 0.385 Figure 4: Linear model: the brake torque is independent of angular velocity Figure 5: Biased exponential model: the braking torque is independent of the square of angular velocity 39 n0⋅N⋅z / 60, that is in case of n0 = 4000 rpm and N = 100 samples, fS = 440000 Hz. The nearest possible setting of the DSO is 500 kHz. The capacity of the sample storage is 50000 samples, so the duration of the measurable signal is 0.1 s. As it is greater than the duration T2 of two crankshaft revolution, the measurement of at least one full engine cycle is feasible using this DSO (T2 = 2⋅(60/n0) = 0.03 s). The instantaneous angular velocity is the reciprocal of the time between two falling edge of the sensor’s signal, so the raw signal had to be smoothed before the zero-crossing detection. Two different smoothing methods were implemented in this case as well: a 10-point moving average smoothing and SAVITZKY-GOLAY smoothing with tenth order polynomial and 31-point length window (for the latter see Fig.7). Because the number of teeth is known, the exact value of the crankshaft angle also with the tooth number can be tracked (see Fig.8). As the angular velocity values are non- equidistant, cubic spline interpolation-based re- sampling was applied as in the previous section. Due to the fast sampling rate, we can track very short-time changes in angular velocity values during 100 ms. However, our task is to estimate the crankshaft’s instantaneous angular velocity. For this purpose the trend of the angular velocity curve estimated at teeth is necessary. The trend was determined on the one hand by using a 10-point moving average filter, on the other with FFT-based filtering (for the latter see Fig.9). In the first case the time of averaging was around 2 ms. Using this value for the time-constant τ of a low-pass filter for the FFT-based computation, for the cut-off frequency of the low-pass filter we get fH = 1 / (2⋅π⋅t) ≈ 80 Hz. Following its definition, the running irregularity was computed using the maximum, minimum and average values of the angular velocity trend shown in Eq.(14). δ = ωmax −ωmin ωaverage = 389.3−379.7 384.7 ≈ 0.02 . (14) The estimation of instantaneous acceleration was computed with differences in case of moving average- based trend. In case of FFT-based trend, the SAVITZKY- GOLAY smoothing derivative algorithm was applied (tenth order polynomial, 31-point length window). Fig.10 illustrates the instantaneous acceleration as a function of time. In this measurement the maximum Figure 6: Setup for measuring the instantaneous torque Figure 8: Inductive sensor’s signal labelled with the crankshaft’s angle and tooth number (angles are in degrees) Figure 7: The raw and smoothed signal of inductive sensor in case of SAVITZKY-GOLAY smoothing Figure 9: Instantaneous angular velocity values and the FFT- based trend Figure 10: Estimation of instantaneous acceleration as smoothing derivative of the trend of instantaneous angular velocity 40 value of angular acceleration was 1113 rad s-2, so for the corresponding instantaneous torque we get: Mmax =Θβmax ≈ 0.004 kg ⋅m2 ⋅1113 rad s2 ≈ 4.5N ⋅m. (15) Nearby the power stroke instant higher value of angular acceleration can be estimated. In that measurement the corresponding instantaneous torque was of 16 N·m (for details see Ref. [6]). The crankshaft’s torque gets to wheel by a drive train with a reduction ratio of 7, so the maximal torque at the wheel is 112 N·m. The wheel’s maximal slipless torque transmission is 40–50 N·m that is during the maximal instantaneous acceleration the wheel will slip on the road, which results in energy loss. This effect ought to be concerned in designing the drive chain of Megameter-IV, e.g. torque-damping torsion clutch should be built in. Conclusions Two methods were elaborated for instantaneous speed and acceleration of an internal combustion engine via on-line measurement. Both the injector-pulse-based and the inductive sensor-based methods gave useful data for developing the next engine, namely the Megameter-IV. The results presented in this paper are the numerical estimation of the engine’s internal friction, the damping constant, the instantaneous angular velocity, the instantaneous angular acceleration and the estimation of the instantaneous torque of the engine during the power stroke. These results provided the foundation for the design and construction of the on-board electronic control unit, which can be seen on Fig.11. Acknowledgement This research is supported by TÁMOP-4.2.2.C- 11/1/KONV-2012-0012: "Smarter Transport" - IT for co-operative transport system and TÁMOP-4.2.2.A- 11/1/KONV-2012-0012: Basic research for the development of hybrid and electric vehicles, The Project is supported by the Hungarian Government and co- financed by the European Social Fund. REFERENCES [1] www.eco-marathon.eu (last accessed: May 20, 2014) [2] BÉLANGER P.R.: Estimation of angular velocity and acceleration from shaft encoder measurements, Proc. IEEE International Conference on Robotics and Automation, Nice, France, 1992, 585–589 [3] STOTSKY A., FORGO A.: Recursive spline interpolation method for real time engine control applications, Control Engineering Practice, 2004, 12, 409–416 [4] CRUZ-PERAGÓN F., PALOMAR J.M., DÍAZ F.A., JIMÉNEZ-ESPADAFOR F.J.: Fast on-line identification of instantaneous mechanical losses in internal combustion engines, Mechanical Systems and Signal Processing, 2010, 24, 267–280 [5] BAGÁNY M.: Internal combustion engines, Kecskemét College, GAMF Faculty, 2011, 192 (in Hungarian) [6] PINTÉR I., BAGÁNY M.: Algorithm development for the Megameter-III telemetry system to estimate the parameters of extremely low fuel consumption internal combustion engine, Proc. Factory Automation, Veszprém, Hungary, 2013, 21–22, 94–99 Figure 11: ECU is under development for motor control and telemetry HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 41–46 (2013) INFLUENCE OF CAN FLATNESS ON HEAT DISSIPATION OF ALUMINIUM ELECTROLYTIC CAPACITOR LÁSZLÓ KOVÁCS1!, LÁSZLÓ GÁL1, AND DÉNES FODOR2 1 Department of Aluminium Electrolytic Capacitor Development (SZ ALU PD), Epcos LLC, Szt. László Str. 6, Szombathely, H-9700 HUNGARY 2 Automotive Mechatronics Department, Institute of Mechanical Engineering, Faculty of Engineering, University of Pannonia, Egyetem Str. 10, Veszprem, H-8200 HUNGARY !E-mail: laszlo.kovacs@epcos.com The lifetime of aluminium electrolytic capacitors highly depends on their core temperature. Heat dissipation in general applications happens by the extended cathode, which is in contact with the inner side of the can. In the case of heat sink applications, the most important heat transfer phenomenon is the heat conduction through the bottom of the aluminium can. The quantity of the dissipated heat is in direct proportion to the size of the heat transfer surface. The more dissipated heat may increase the lifetime of the capacitor. Therefore, the flatness value of the can bottom is critical. This paper presents a flatness measurement method, which can successfully replace the equipment for a more complex and more expensive 3D measurement. It discusses an implementation of a measurement environment, where data acquisition and visualization are automated by a LabVIEW-based software. In addition, this study deals briefly with the influence of production processes on the flatness value of the capacitor produced by leading manufacturers. Keywords: aluminium electrolytic capacitor, measurement automation, heat dissipation, graphical programming environment, flatness measurement Introduction The aluminium electrolytic capacitor is the most commonly applied capacitor type, due to the fact that the capacitance and voltage range of the components is wide. The rated voltages change from 5 V to 550 V, while the capacitance 1 µF to 3 F. These are applied in many fields of industry, such as energetics, power electronics, automotive application, etc. and used for energy storage, smoothing and filtering function. The lifetime of the capacitor highly depends on its core temperature. Each decrease of 10 ºC doubles the lifetime [1]. Operating temperature is determined by the ambient temperature, by the applied ripple current, by the used voltage, and by the equivalent serial resistance (ESR). There are two ways for the reduction of core temperature. The first one is the extended cathode foil and the second one is the cooling of the capacitor can. Cooling in general applications is realized by applying a heat sink at the bottom of a capacitor can. This paper presents the basic construction of capacitors and the most important production steps, introduces the heat conduction between the can and the environment and describes the effect of the flatness of the can to heat conductivity. It also presents the entire measurement environment (measurement station, data acquisition and evaluation software), the database behind the measurement system that stores the results. At last, it discusses the effect of production processes to the measured flatness, and presents the flatness of capacitors produced by leading capacitor manufacturer companies. Structure and construction of aluminium electrolytic capacitor The winding of an aluminium electrolytic capacitor contains two foils and papers [1]. These are rolled together tightly into a winding. The material of the anode, positive foil is aluminium with purity higher than 99.9%. The foil has been etched [2] to increase the effective surface area (and thus the capacitance of the capacitor). As a result, the effective surface area becomes typically 20–40 times larger than the plain area of the foil. On the etched surface of the foil an aluminiumoxide layer [3] has been generated electrochemically. The forming voltage of the anodized aluminium foil [4] is 30–60% higher than the rated voltage of the capacitor. The material of the cathode foil is also aluminium and itself has a thin oxide film (the forming voltage is only a few volts regardless of the rated voltage). It is typically etched to slightly increase the surface area. The anode and cathode foils are connected to aluminium tabs, which are coming out from the winding and are riveted to the aluminium terminals of the cover disk. The tab foils are not etched but are provided with an oxide layer made by 42 electrochemical oxidization. Before being housed in a suitable container, the complete winding is impregnated with electrolyte. After housing the edges of the can are curled down. Before being sleeved and packed, capacitors are first aged. The purpose of this stage is to repair any damage in the oxide layer and thus to reduce the leakage current to very low levels. During manufacturing there are two processes (curling and aging processes), which apply mechanical stress to the casing of capacitors. The effect of flatness of aluminium cans on heat conductivity In common applications with usage of extended cathode (the cathode is wider than the anode and exposed to the bottom of the winding) or different kinds of capacitor cooling (for example air or water cooling) are applied to reduce the core temperature that effects its lifetime. Earlier investigations showed that there are three types of heat transfer phenomenon between capacitor and the environment [5]: (i) heat conductivity (mostly at the bottom), (ii) radiation (mostly at the side), and (iii) convection. This paper only deals with heat transfer from the bottom of the capacitor to the environment. From the general theory of thermodynamics the heat conduction can be described by Eq.(1). P = kA dT dx , (1) where P represents power of the heat flow (W); A is the size of the tangential surface (m2); dT/dx is the temperature gradient (K m-1); and k stands for heat conductivity of the material (W Km-1). In our case the linear approximation is acceptable. As already described, the rate of the heat flow depends on the size of the tangential plane. Definition of the flatness was adapted from Ref. [6]. It is given as the distance between the two closest tangential planes of the bottom of the aluminium electrolytic capacitor can. In the case of capacitors with higher flatness the tangential surface of the aluminium can is lower. The distance between the bottom of the capacitor and the heat sink is filled by air. It is well known that the heat conduction coefficient of the air is much lower than the aluminium. The aspect ratio between the two constants is in the order of magnitude ~104. Therefore the flatness of the bottom in the case of regular heat sink applications is an important value. Theoretical background of used measurement method There are several different methods for measuring the flatness of a surface [4]: 1. The entire surface is compared with a known reference surface (holistic methods); 2. Points on the surface are compared to a reference plane; 3. Straightness of the lines in the surface is measured. The holistic methods (e.g. holographic and interferometric methods) are not suitable for large surfaces with more than 90 mm diameter of screw terminal capacitor [7]. The third method is too complex and time-consuming for applying in serial production. Therefore, a point-to-point method (second type) was chosen for this investigation. Of course, the expensive optical methods may be more precise (in the order of 100 nm deviations) [8], but in our case this high precision is unnecessary. On the other hand, the point- to-point method is very fast, can be used easier and cost-efficient. For point-by-point methods, the grid size is a critical issue, which determines the lowest limit for observable surface irregularities. In our case, approximately 60 cm2 surface area needs to be characterized. The height difference was measured between N points within the bottom. Trials were made to find the applicable grid. Finally, a 9-point-grid was chosen with circular symmetry, because the measurement time is short enough and the accuracy of the measurement is acceptable. Least mean squares (LMS) method was used for plane fitting, because the randomly distributed measurement errors have the smallest influence to the results [9]. Eq.(2) was used for the LMS: xi 2 i=1 N ∑ xiyi i=1 N ∑ xi i=1 N ∑ xiyi i=1 N ∑ yi 2 i=1 N ∑ yi i=1 N ∑ xi i=1 N ∑ yi i=1 N ∑ 1 i=1 N ∑ " # $ $ $ $ $ $ $ $ % & ' ' ' ' ' ' ' ' × a b c " # $ $ $ % & ' ' ' = zixi i=1 N ∑ ziyi i=1 N ∑ zi i=1 N ∑ " # $ $ $ $ $ $ $ $ % & ' ' ' ' ' ' ' ' , (2) where N is the number of the grid points; xi, yi are coordinates within the measured surface [mm]; zi corresponds to the measured height compared to z1 of the ith point (mm), a, b, c are constants for the fitting. The fitted plane is described by Eq.(3). z(x, y) = ax+by+c (3) By solving the system of linear equations, for example with Gauss-elimination the fitted plane is given. The Gauss-elimination is not necessary if the grid points are centrally symmetric. In this case, the input matrix is diagonal, so the solution of the equation system is simple. The distance between the individual points and the fitted plane was calculated. The sign of this distance describes the location of the point respect to the fitted plane. The flatness is the distance between the farthest points above and under the fitted plane. The Fig.1 shows an example of measured points and the fitted plane. 43 Implemented measurement environment The base of the hardware was a measuring pad made from stainless steel, which has an extremely smooth surface. By the LMS method, the error that comes from the tilt of the steel pad was avoided. To stay clear of the effect of the temperature gradient, the measuring equipment was kept under controlled environmental conditions. The data acquisition programme for the flatness measurement was written in LabVIEW with contains two major parts. The first one collects the data from the controlled hardware (height meter) and stores them into a local database in MS Excel format. This software communicates through an RS-232 port with the measuring equipment. The user sets the actual position, pushes a button and the module registers the results. The graphical user interface developed is shown in Fig.2 with Hungarian annotations. The software displays the number of actual measurement points and controls the entire process. For example, if the flatness value exceeds a predetermined limit, the software sends a warning to the operator. The second module can extract the results from the database of an examined capacitor(s) or the user can import the measured values from Excel. If the results are in the input table, the evaluation module can calculate the flatness of the capacitor can. The graphical user interface of the second module can be seen in Fig.3. The implemented automation environment fits into the measurement automation system (MAS) [10] being used at the Aluminium Electrolytic Capacitor Development Department of Epcos LLC. The system uses two frameworks: a LabVIEW based framework for data acquisition and an ASP.NET based framework for data management. The measurement part contains automated data acquisition measurements, which are connected to the electrolyte (conductivity, pH, viscosity, etc.) and different kind of capacitor tests. (Lifetime, surge voltage and storage test, etc.). The experiments on the electrolyte are controlled by an NI-PXI, which is connected to the database (Fig.4). PXI stand for PCI eXtensions for Instrumentation. These platforms are used as a basis for building of electronic test equipment, automation system, etc. PXI chassis can handle many modules for example plug-in data acquisition (DAQ) cards, communication cards like RS-232, and different analogue and digital I/O boards. The capacitor measurements are controlled by computers and not by the NI-PXI, because the current implementation measures the low and high voltages separately. The ASP.NET part of the system was developed for data management and evaluation. This is a software module that contains useful tools that facilitate data handling. Data management module simplifies the registration of the constructive properties of the capacitor (like anode foil, cathode foil, type of can, cover disk, etc.) and helps the data storage of applied Figure 1: Example of the measured points and the fitted plane Figure 3: Graphical user interface of the evaluation software Figure 2: Graphical user interface of the data acquisition software Figure 4: Architecture of the electrolyte measurements 44 voltage and current, ambient temperature, etc. The evaluation part supports the evaluation process by generating a standardized report. The user can tailor the reports according to its needs by the Report Generation tool. The desired data appear in a representative way and even the trends of the parameters can be shown. The goal of the MAS is to automate the previously manual measurements and eliminate paper-based registration. There are many advantages like making the measurements more precise, more reliable and fault tolerant, running multiple measurements in parallel, which all contribute to speed up the research and development of new component and devices. Change of flatness value during the manufacturing The flatness of the capacitor can is measured at the incoming inspection and at the final measurement of the production. The initial value of the flatness of the capacitor is determined by the incoming aluminium cans. These parts are produced by cold extrusion of aluminium slugs, so the irregularities of the cans is in the order of magnitude ~10 µm, which is negligible for our purposes. The Fig.5 shows the flatness values measured at the incoming inspection of the capacitors. It can be seen that the mean of the measured flatness values is 0.09 mm. As mentioned above, the flatness is mainly affected by two procedures: curling and aging. During curling the can is closed hermetically by curling back the edges of the aluminium can. This curled edge sinks into the rubber ring of the cover disk. As a consequence, the cover disk presses down the winding. It can affect the flatness value, since for example concave shaped aluminium can turn into a convex shape. During aging, the flaws of the anode foil are repaired by applying voltage to the capacitor and placing it in an oven. As a side effect gas is generated. The pressure of the generated gas can change the geometry of the can. From the distribution of the flatness values the effect of Epcos production steps (coloured curves in Fig.6) can be estimated. The results show that the different production steps do not influence significantly the bottom flatness of the capacitor at the Epcos. Experimental results regarding to core temperature measurements Two experiments were completed for investigating the importance of the flatness value of the capacitor can. Sample capacitors were produced with different bottom flatness values. The aluminium cans were made to be convex. These were ordered directly from the supplier and the used parts were sorted out. Regular heat sinks were assembled to the bottom of the parts. The tangential planes of the heat sinks were flat (~2 µm) due to their grinded surface. The heat sinks were fixed to the bottom of the samples with the bottom screw of the capacitors. The mechanical stability of the fixation was made with torque wrench. Also guaranteed the same fixation level in each case. During the tests, forced air- cooling was not used. Heat conductive wires were applied to the core of the samples, which allowed us to measure the core temperature of the capacitors directly. The core temperature measurement system is shown in Fig.7. In case of the first measurement, sample capacitors with very similar equivalent serial resistance (ESR) and leakage current were selected. These parts were taken to 85 °C oven and rated voltage (450 V) was applied to them for 100 h. This preparation is necessary to avoid the different behaviour of the leakage current of the capacitors (isolated from each other to avoid the heat transfer). During the experiment a sinusoidal current Figure 5: The measurement results from Incoming Inspection Department Figure 6: The effect of the production processes to the flatness of capacitor Table 1: Results of the ripple current measurement Average Flatness (mm) Average of Temperature at equilibrium state (°C) 0.102 90.6 0.259 90.9 0.509 91.5 45 was applied at 50 Hz. The current and the temperature (the upper category temperature) of the test are specified by the data sheet of the capacitor. In this case the examined capacitor was B43586A5278Q (C = 2700 µF, UR = 450 V, I (85 °C, 100 Hz) = 12 A), the applied current (multiplying with the frequency factor) was 10.8 A. The selected samples were placed into an oven at 85 °C and a rated voltage of 450 V was applied. The core temperature was measured during the test. The equilibrium temperature was measured and compared. Test was started with the remaining three prepared samples. In case of the second measurement the behaviour of the capacitors was tested against transient heat. The parts were heated up to 85 °C and removed from the oven after 10 hours. The core temperature of the capacitors was measured during cool down time. The recorded core temperature curves (Fig.8) show exponential decay. The time-constant of the decay depends on the flatness value of the bottom. As a comparison, the flatness values for screw terminals capacitors from leading manufacturers were inspected. Table 2 summarizes the flatness results by leading manufacturers. Conclusion It is possible to measure the flatness of a capacitor can with the presented and implemented measurement environment. The expensive 3D measurement equipment can be replaced with this cost-effective method. The measurement and evaluation software was implemented in LabVIEW programming environment, which perfectly fits into the previously presented measurement automation system, which is currently being used at the aluminium electrolytic capacitor product development department of Epcos. The faults of previously manual measurements have been lowered with the help of automation and data acquisition software. This study discussed the basic construction of a capacitor and highlighted the critical construction steps from the aspect of flatness. In addition, it reveals the flatness values of other capacitor manufacturers’ products. Acknowledgments The authors gratefully acknowledge the support of the Hungarian State and European Union under the "Green Energy – Cooperation of the higher education sector for the development of green economy in the area of energetics” project (TÁMOP-4.1.1.C-12/1/KONV- 2012-0017), and the support of the TÁMOP-4.2.2.A- 11/1/KONV-2012-0071. REFERENCES [1] TDK Corporation, Aluminium Electrolytic Capacitors, 2013 [2] PARK D., KIM H.: Electrochemical Etching of Aluminium through Porous Alumina, Analytical sciences, 2001, 17(supplement), a73–a76 [3] CHURCH H.F.: The dielectric properties of anodic aluminium oxide films, Proc. IEE – Part B: Electronic and Communication Engineering, 1962, 109, 399–403 Figure 8: The behaviour of core temperature in case of different bottom flatness values Figure 7: Core temperature measurement of a capacitor Table 2: Flatness results for products from leading manufacturers Company Capacitor Series Flatness (mm) Epcos B43586 0.095±??? Hitachi FXR 0.158±??? Kemet ALS 0.225±??? Nichicon NT 0.330 46 [4] CHEN C., MASS W., HUTCHINS G.: Anodization of aluminium electrolyte capacitor foil, U.S. Patent 4,481,084, 1984 [5] PARLER S.G. JR.: Thermal Modelling of Aluminium electrolytic Capacitor, IEEE Industry Applications Society Conference, 1999 [6] ISO/TS 12781–1: Geometrical Product Specifications (GPS) – Flatness – Part 1: Vocabulary and parameters of flatness, 2003 [7] DE GROOT P.J.: Grating interferometer for flatness testing, Optics Letters, 1996, 21(3), 228–230 [8] MEIJER J.: Accuracy of Surface Plate Measurements – General Purpose Software for Flatness Measurement, Annals of the CIRP, 1990, 39, 545–548 [9] HAITJEMA H., MEIJER J.: Evaluation of surface plate flatness measurements, European Journal of Mechanical Engineering, 1993, 38(4), 165–172 [10] KOVÁCS L., FODOR D., KLUG O., ENISZ K.: Measurement Automation System for Aluminium Electrolyte Capacitor Development, IEEE Instrumentation & Measurement Magazine, 2013, 16(3), 38–43 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 51-57 (2013) QUASI-POLYNOMIAL CONTROL OF A SYNCHRONOUS GENERATOR ATTILA MAGYARB1 AND ATTILA FODOR1 1Department of Electrical Engineering and Information Systems, University of Pannonia, Egyetem u. 10., H8200 Veszprém, HUNGARY BE-mail: magyar.attila@virt.uni-pannon.hu A simple dynamic model of permanent magnet synchronous generator, that is used for electrical energy generation is investigated in this work using a nonlinear technique based on the quasi-polynomial representation of the dynamical model. It is well known that dynamical systems with smooth nonlinearities can be embedded in a quasi-polynomial model. Quasi- polynomial systems are good candidates for a general nonlinear system representation since their global stability analysis is equivalent to the feasibility of a LMI. Moreover, the stabilizing quasi-polynomial state feedback controller design problem is equivalent to the feasibility of a bilinear matrix inequality. The classical stabilizing state feedback problem for quasi- polynomial systems has been extended in this work with the ability of tracking time-dependent reference signals. It is shown, that the stabilizing quasi-polynomial servo controller design is equivalent to a bilinear matrix inequality. The results are applied to the model of a synchronous generator. Keywords: quasi-polynomial systems, Lotka-Volterra systems, stability analysis, state feedback control, synchronous generator, wind turbine Introduction Electrical power systems should operate in an economic way with minimum possible operating cost under nor- mal operating conditions. To ensure this, a preventive controller for power systems has been presented in [1]. It encompasses many types of control actions, including generation rescheduling, load curtailment and network switching reactive compensation. From the viewpoint of the power grid, the electric power generation can be characterized by the operation of the electrical generators, the subject of our study. These power plants should not only be able to follow the time- varying active power demand of the consumers and the central dispatch center, but also keep the quality indica- tors (frequency, waveform, total harmonic distortion) of the grid on the expected level. This can be achieved by applying proper control methods based on dynamic mod- els of plant (see e.g. [2, 3]) and the involved generators. Because of the specialties and great practical impor- tance of the synchronous generators in power plants, their modeling for control purposes is also well investigated in the literature. Besides of the basic textbooks (see e.g. [4]) that describe the modeling, specialized papers are also available that use the developed models for the design of various controllers [5]. A wind turbine driving permanent magnet syn- chronous generator is proposed in [6] with current con- trolled voltage source inverter, which is the best choice when the output power is small. The current control of the voltage source inverter has bidirectional active and reactive power control ability, avoiding the intricacy of the controller designing. The class of quasi-polynomial (QP) systems plays an important role in the theory of nonlinear dynamical sys- tems because nonlinear systems with smooth nonlineari- ties can be transformed into a QP form [7]. This means, that any applicable method for QP systems can be re- garded as a general technique for nonlinear systems. Previous work in the field of QP systems include [8], which proves that the global stability analysis of QP sys- tems is equivalent to the feasibility of a linear matrix in- equality (LMI). It has been shown in [9] that the globally stabilizing state feedback design for QP systems is equiv- alent to a bilinear matrix inequality (BMI). Although the solution of a BMI is an NP-hard problem, an itera- tive LMI algorithm could be used. A summary of linear and bilinear matrix inequalities and the available software tools for solving them can be found in [10]. Another con- trol synthesis algorithm for polynomial systems is pre- sented in [11]. The goal of this paper is to formulate the servo con- troller design problem for QP systems based on the re- sults presented in [9] and to design a servo controller for a synchronous generator model using the QP controller synthesis methodology that keeps the active power at the desired level. 52 Basic notions In what follows, the basic modeling assumptions and def- initions to be used in the sequel are summarized briefly. Nonlinear model of a synchronous generator The modeling procedure of the synchronous generator is mainly based on [12] and [13], therefore, only the result- ing model is presented here. The model is based on the following simplification assumptions: • a symmetrical tri-phase stator winding system is as- sumed, • one field winding is considered to be in the machine, • all of the windings are magnetically coupled, • the flux linkage of the winding is a function of rotor position, • the copper loss and the slots in the machine can be neglected, • the spatial distribution of the stator fluxes and aper- tures wave are considered to be sinusoidal, and • stator and rotor permeability are assumed to be infi- nite. It is also assumed that all the losses due to wiring, satu- ration, and slots can be neglected. The four windings (three stators and one rotor) are magnetically coupled. Since the magnetic coupling be- tween the windings is a function of the rotor position, the flux linkage of the windings is also a function of the rotor position. The actual terminal voltage v of the windings can be written in the form v = ± JX j=1 (r j i j )± JX j=1 ✓ d' j dt ◆ , (1) where i j are the currents, r j are the winding resistances, and ' j are the flux linkages. The positive directions of the stator currents point out of the synchronous generator terminals. Thereafter, the two stator electromagnetic fields, both traveling at rotor speed, can be identified by decompos- ing each stator phase current under steady state into two components, one in phase with the electromagnetic field and another phase shifted by 90�. With the above, one can construct an air-gap field with its maximal aligned to the rotor poles (d axis), while the other is aligned to the q axis (between poles). This method is called the Park’s transformation. As a result of the derivation in [12] the vector voltage equation is as follows v dFq = �Ri dFq �L d dt i dFq , (2) with v dFq = ⇥ v d �v F v q ⇤ T i dFq = ⇥ i d i F i q ⇤ T , (3) Figure 1: The equivalent circuit of the synchronous generator. where v d , v q and i d , i q are the direct and the quadrature components of the stator voltage and current of the syn- chronous generator, while v F and i F are the exciter volt- age and current of the synchronous generator. Further- more, R and L are the following matrices (see Fig.1) RRS! = 2 4 r 0 !L q 0 r F 0 �!L d �! kM F r 3 5 L = 2 4 L d kM F 0 kM F L F 0 0 0 L q 3 5 , (4) where r is the stator resistance, r F is the exciter resis- tance of the Synchronous Generator, L d , and L q are the direct and the quadrature part of the stator and rotor in- ductance, ! is the angular velocity, and M F is linkage inductances. The state-space model for the currents is ob- tained by expressing d dtidFq from Eq. (2), i.e. d dt i dFq = �L�1 R RS! i dFq �L�1 v dFq . (5) The purely electrical model Eq. (5) has to be extended with the equation of rotational motion (Eq. (6)) that gives the mechanical sub-dynamics, that is d! dt = �L d i q 3⌧ j i d + �kM F i q 3⌧ j i F + + L q i d 3⌧ j i q + �D ⌧ j ! + Tmech ⌧ j . (6) 2 664 i̇d ˙ iF i̇q !̇ 3 775 = 2 6666666664 �LF H kMF H 0 0 kMF H �Ld H 0 0 0 0 1 Lq 0 0 0 0 1 ⌧j 3 7777777775 2 64 �vd vF �vq Tmech 3 75+ 53 2 66666666664 � rLF H kMF rF H �!LF Lq H 0 rkMF H � rF Ld H !kMF Lq H 0 �!Ld Lq !kMF Lq r Lq 0 �Ld iq 3⌧j �k MF iq 3⌧j Lq id 3⌧j �D ⌧j 3 77777777775 2 64 id iF iq ! 3 75 (7) where H = k2M2 F � L d L F 2 664 i̇d ˙ iF i̇q !̇ 3 775 = 2 64 �1.7100 0.5893 0 0 0.5893 �6.6918 0 0 0 0 �1.7090 0 0 0 0 0.0006 3 75 2 64 �vd vF �vq Tmech 3 75+ 2 64 �0.0019 0.0004 �3.4883! 0 0.0006 �0.0042 1.2022! 0 3.5888! 2.6489! �0.0019 0 �0.0004 iq �0.0003 iq 0.0004 id �0.0011 3 75 2 64 id iF iq ! 3 75 (8) Altogether, the state equations Eqs. (5) and (6) have four state variables: i d , i F , i q , and !. The manipulated input vector of the generator is u = ⇥ v F T mech ⇤ T , while the disturbance input vec- tor is d = ⇥ v d d q ⇤ T . Realize that the state equations Eqs. (5) and (6) are bilinear in the state variables because matrix RRS! depends linearly on !. The obtained model is nonlinear and it has four state variables: i d , i F , i q and ! [14]. Output equations of the model: The output active power equation can be written in the following form: p out = v d i d + v q i q (9) and the reactive power is q out = v d i q � v q i d . (10) Eqs. (9) and (10) are the output equations of the genera- tor’s state space model. It is important to note, that these equations are bi-linear in the state and input variables. Note, that although only the active power is to be con- trolled in this case, as a possible future extension of the work, reactive power can also be controlled in order to follow an external reference signal. Quasi-polynomial representation of nonlinear systems Let us denote the element of an arbitrary matrix W with row index i and column index j by W ij . Quasi- polynomial models are systems of ordinary differential equations (ODE) of the following form ẏ i = y i 0 @L i + mX j=1 A ij nY k=1 y Bjk k 1 A , i = 1, . . . , n. (11) where y 2 int(Rn +), A 2 Rn⇥m, B 2 Rm⇥n, L i 2 R, i = 1, . . . , n. Furthermore, L = [L1 . . . L n ]T . Let us denote the equilibrium point of interest of Eq. (11) as y⇤ = [y⇤1 y⇤2 . . . y⇤ n ]T . Without the loss of generality we can assume that rank(B) = n and m � n (see [15]). Lotka-Volterra models The above family of models is split into classes of equivalence [16] according to the values of the products M = BA and N = BL. The LOTKA-VOLTERRA form gives the representative elements of these classes of equivalence. If rank(B) = n, then the set of ODEs in Eq. (11) can be embedded into the following m- dimensional set of equations, the so called LOTKA- VOLTERRA model: ż j = z j N j + mX i=1 M ji z i ! , j = 1, . . . ,m (12) where M = BA, N = BL, and each z j represents a so called quasi-monomial: z j = nY k=1 y Bjk k , j = 1, . . . ,m. (13) Stability analysis using linear matrix inequalities Henceforth it is assumed that y⇤ is a positive equilib- rium point, i.e. y⇤ 2 int(Rn +) in the QP case and simi- larly z⇤ 2 int(Rm + ) is a positive equilibrium point in the LOTKA-VOLTERRA case. For LV systems there is a well known candidate LYAPUNOV function family ([8, 17]), which is in the form: V (z) = mX i=1 c i ✓ z i � z⇤ i � z⇤ i ln z i z⇤ i ◆ , (14) c i > 0, i = 1 . . .m, where z⇤ = ⇥ z⇤1 . . . z⇤ m ⇤ T is the equilibrium point corresponding to the equilibrium y⇤ of the original QP system (Eq. (11)). The time derivative of the LYAPUNOV function Eq. (14) is: V̇ (z) = 1 2 (z � z⇤)(C M +MTC)(z � z⇤), (15) where C = diag(c1, . . . , cm) and M is the invariant characterizing the LOTKA-VOLTERRA form (Eq. (12)). Therefore, the non-increasing nature of the LYAPUNOV function is equivalent to a feasibility problem over the following set of LMI constraints (see [18] or [19]): C M +MTC  0 C > 0 (16) where the unknown matrix is C, which is diagonal and contains the coefficients of Eq. (14). 54 The derivation of global stability analysis for nonau- tonomous QP systems from the autonomous case is straightforward. The LYAPUNOV function (Eq. (15)) also depends on the equilibrium value of the input (u⇤) and has the form ˙̃V (z) = 1 2 (z � z⇤)(C M̃ + M̃ T C)(z � z⇤), (17) where M̃ depends on the coefficient matrices of the input-affine LOTKA-VOLTERRA model (Eq. (21)): M̃ = M0 + pX l=1 Ml u ⇤ l . The corresponding LMI feasibility problem to be solved in order to check global asymptotic stability is C M̃ + M̃ T C  0 C > 0. (18) Input-affine QP system models An input-affine nonlinear system model with state vector y, input vector u and output vector ⌘ ẏ = f(y)+ P p i=1 gi(y)ui ⌘ = h(y) (19) is in QP-form if all of the functions f , gi and h are in QP-form. Then the general form of the state equation of an input-affine QP system model with p-inputs is: ẏ i = y i 0 @L0i + mX j=1 A0ij nY k=1 y Bjk k 1 A+ (20) + pX l=1 y i 0 @L li + mX j=1 A lij nY k=1 y Bjk k 1 Au l where i = 1, . . . , n, A0,Al 2 Rn⇥m, B 2 Rm⇥n, L0,Ll 2 Rn, l = 1, . . . , p. The corresponding input-affine LOTKA-VOLTERRA model is in the form ż j = z j N0j + mX k=1 M0jkzk ! + (21) + pX l=1 z j N lj + mX k=1 M ljk z k ! u l where j = 1, . . . ,m, M0,Ml 2 Rm⇥m, N0,Nl 2 Rm, l = 1, . . . , p, and the parameters can be obtained from the input-affine QP system’s ones in the following way M0 = BA0 N0 = BL0 Ml = BAl Nl = BLl l = 1, . . . , p. (22) The controller design problem for QP systems Globally stabilizing QP state feedback design problem for QP systems can be formulated as follows (for a more detailed desription, see [9]). Consider arbitrary QP inputs in the form: u l = rX i=1 k il q̂ i , l = 1 . . . , p, (23) where q̂ i = q̂ i (y1, . . . , yn), i = 1, ..., r are arbi- trary quasi-monomial functions of the state variables of Eq. (20) and k il is the constant gain of the quasi- monomial function q̂ i in the l-th input u l . The closed loop system will also be a QP system with matrices  = A0 + pX l=1 rX i=1 k il Ail, B̂, L̂ = L0 + pX l=1 rX i=1 k il Lil. (24) Note that the number of quasi-monomials in the closed- loop system (i.e. the dimension of the matrices) together with the matrix B̂ may significantly change depending on the choice of the feedback structure, i.e. on the quasi- monomial functions q̂ i . Furthermore, the closed loop LV coefficient matrix M̂ can also be expressed in the form M̂ = B̂  = M0 + pX l=1 rX i=1 k il Mil. (25) Then the global stability analysis of the closed loop sys- tem with unknown feedback gains k il leads to the follow- ing BMI M̂ T C +CM̂ = M0 TC +CM0+ (26) pX l=1 rX i=1 k il ⇣ Mil TC +CMil ⌘  0. The variables of the BMI are the p⇥ r k il feedback gain parameters and the c j , j = 1, ..,m parameters of the LYAPUNOV function. If the BMI above is feasible, there exists a globally stabilizing feedback with the selected structure. Note that (marginal) infeasibility of the BMI (Eq. (26)) means only that the closed loop system is not proven to be globally asymptotically stable. However, the solution k il may still guarantee local stability, which is enough in several cases. Controller design using bilinear matrix inequalities A BMI is a diagonal block composed of q matrix inequal- ities of the following form Gi 0 + pX k=1 x k Gi k + pX k=1 pX j=1 x k x j Ki kj  0, i = 1, . . . , q (27) 55 where x 2 Rp is the decision variable to be determined and Gi k , k = 0, . . . , p, i = 1, . . . , q and Ki kj , k, j = 1, . . . , p, i = 1, . . . , q are symmetric, quadratic matrices. The main properties of BMIs are that they are non- convex in x (which makes their solution numerically much more complicated than that of linear matrix in- equalities), and their solution is NP-hard [10], so the size of the tractable problems is limited. However, there exist practically applicable and effective algorithms for BMI solution [20], [21], or [22]. In Matlab environment the TomLab/PENBMI solver [23] can be used effectively to solve BMIs. Quasi-polynomial servo control As in the linear case, the problem statement of the servo, or reference tracking control is as follows. Consider a nonlinear system in the form of Eq. (20) and an external reference signal r which is to be followed by the system output ⌘. It is possible to define the tracking error signal (Eq. (28)) whose time derivative gives the tracking error dynamics which should be stabilized together with the system (Eq. (20)). z(t) = Z t t0 r(⌧)� ⌘(⌧)d⌧. (28) It is easy to see that the differential equation of the track- ing error has the form (Eq. (29)). ż(t) = r(t)� ⌘(t), (29) u(t) = �K R Z t 0 r(⌧)� ⌘(⌧)d⌧ �K y(t). If the output equation is also in QP form, then the ex- tended closed loop QP system can be written up in LOTKA-VOLTERRA form similar to (Eq. (25)) and the BMI for the globally stabilizing controller design can be formulated. Quasi-polynomial control of the synchronous generator Quasi-polynomial form of the synchronous generator The bilinear nature of the state equations (Eq. (8)) enables us to directly apply the QP mechanism without QP em- bedding [16]. The system has the following set of quasi- monomials: ⇢ 1 i d , i F i d , i q ! i d , 1 i F , i d i F , i q ! i F , 1 i q , i d ! i q , i F ! i q , 1 ! , i F i q ! � The QP coefficient matrices of the input-affine system (Eq. (20)) are A0 = 2 66666666666666664 �2.323 0 0 0 0.0004 0 0 0 0.5893 0 0 0 �3.4883 0 0 0 0 0.8 0 0 0 0.0006 0 0 0 �6.6918 0 0 0 1.2022 0 0 0 0 �0.7861 0 0 0 3.5888 0 0 0 2.6489 0 3 77777777777777775 T A1 = 2 66666666666666664 0.5893 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 �6.6918 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 77777777777777775 T A2 = 2 66666666666666664 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0006 3 77777777777777775 T Feedback structure There is a degree of freedom in the selection of the stabi- lizing feedback structure. A wise choice of the feedback structure does not increase the number of monomials of the closed loop system. This way the size of the BMI to be solved remains tractable. In our case, a linear full state feedback is applied, i.e. the feedback law is in the form v F = k r1 R t 0 r(⌧)� p out (⌧)d⌧+ +k1 id + k2 iF + k3 iq + k4 ! T mech = k r2 R t 0 r(⌧)� p out (⌧)d⌧+ +k5 id + k6 iF + k7 iq + k8 !. (30) Controller design and verification via simulation Using the feedback law (Eq. (30)), the closed loop system is also in QP form with 18 quasi-monomials: ⇢ 1 i d , i F i d , i q i d , ! i d , i q ! i d , 1 i F , i d i F , i q i F , ! i F , i q ! i F , 1 i q , i d ! i q , 56 0 10 20 30 40 50 60 70 80 90 100 0 0.5 1 1.5 time pout(t) r(t) Figure 2: The controlled output p out of the system (solid line) together with the reference input r (dashed line). i F ! i q , i d ! , i F ! , i q ! , i d i q ! , i F i q ! � . Due to the lack of space, the 19⇥19 (the extra dimension comes from the tracking error dynamics) matrices of the BMI and the closed loop LOTKA-VOLTERRA system are not listed here. The BMI (Eq. (26)) for the globally stabilizing feed- back design case suffers from rank deficiency and the available BMI solvers stop with marginal infeasibility, which means that global stability cannot be proven using the feedback law (Eq. (30)). However, simulation results indicate, that the system is locally stable with the con- troller gain parameters yielded by the globally stabilizing BMI (Eq. (26)) formulated for the closed loop generator model extended with the tracking error dynamics. k r1 = 1.0000 k r2 = �0.0002 k1 = �0.2932 k2 = �1.6880 k3 = �1.9609 k4 = 0.2259 k5 = 0.0001 k6 = 0.0000 k7 = �0.0001 k8 = 0.1450 The simulated behavior of the controlled generator can be seen in Fig. 2 where the controlled output of the system (i.e. p out ) is shown together with the reference input (r). Both of them are dimensionless. The tracking properties of the system are acceptable, moreover, its disturbance rejection is also good (at time 40, a step-like change has been applied to the disturbance input v d ). The simulation has been performed in Mat- lab/Simulink environment [24]. Conclusions A novel servo control design technique based on the QP representation has been formulated in this work. As an example, an active power tracking controller has been de- signed for a synchronous generator model. The tracking properties of the closed loop system are satisfactory. Further work includes the extension of the method for vector reference signals and use the degree of freedom ly- ing in the BMI problem for formulating a robust/optimal controller design problem. Another direction of future re- search is to apply graph theoretical techniques for con- troller structure selection that applies the underlying con- nections between QP systems and chemical reaction net- works [25]. Acknowledgments This research was supported by the European Union and the State of Hungary, co-financed by the European Social Fund in the framework of TÁMOP 4.2.4. A/2-11-1-2012- 0001 ’National Excellence Program’. REFERENCES [1] VERMA, K., NIAZI, K. A coherency based gener- ator rescheduling for preventive control of transient stability in power systems, International Journal of Electrical Power & Energy Systems, 2013, 45, 10– 18 [2] BANAVAR, R., DESHPANDE, U. Robust controller design for a nuclear power plant using H-infinity optimization, in Proc. of the 35th Conference on Decision and Control, Kobe, Japan [3] NA, M., JUNG, D., SHIN, S., JANG, J., LEE, K., LEE, Y.: A model predictive controller for loadfol- lowing operation of PWR reactors, IEEE Transac- tions on Nuclear Science, 2005, 52,1009–1020 [4] ANDERSON, P., FOUAD, A. Power-Systems- Control and Stability (The IOWA State University Press, Ames Iowa, USA), 1977 [5] DE MELLO, F., CONCORDIA, C. Concepts of syn- chronous machine stability as affected by excitation control, IEEE Transactions on Power Apparatus and Systems, 1969, PAS-88(4), 316–329 [6] MUYEEN, S., AL-DURRA, A., TAMURA, J. Vari- able speed wind turbine generator system with cur- rent controlled voltage source inverter, Energy Con- vers Manage, 2011, 52 [7] HERNÁNDEZ-BERMEJO, B., FAIRÉN, V. Nonpoly- nomial vector fields under the Lotka-Volterra nor- mal form, Physics Letters A, 1995, 206, 31–37 [8] FIGUEIREDO, A., GLERIA, I. M., ROCHA, T. M. Boundedness of solutions and Lyapunov functions in quasi-polynomial systems, Physics Letters A, 2000, 268, 335–341 [9] MAGYAR, A., SZEDERKÉNYI, G., HANGOS, K. M. Globally stabilizing feedback control of pro- cess systems in generalized Lotka-Volterra form, Journal of Process Control, 2008, 18(1), 80–91 [10] VANANTWERP, J., BRAATZ, R. A tutorial on linear and bilinear matrix inequalities, Journal of Process Control, 2000, 10, 363–385 [11] TONG, C.-F. , ZHANG, H., SUN, Y.-X. Con- troller Design for Polynomial Nonlinear Systems 57 with Affine Uncertain Parameters, Acta Automatica Sinica, 2007, 33(12), 1321–1325 [12] FODOR, A., MAGYAR, A., HANGOS, K. M. Control-Oriented Modeling of the Energy- Production of a Synchronous Generator in a Nuclear Power Plant, Energy, 2012, 39, 135–145 [13] FODOR, A., MAGYAR, A., HANGOS, K. M. Pa- rameter Sensitivity Analysis of a Synchronous Gen- erator, Hungarian Journal of Industrial Chemistry, 2010, 38(1), 21–26 [14] FODOR, A., MAGYAR, A., HANGOS, K. Dynamic modeling and model analysis of a large indus- trial synchronous generator, in Applied Electronics, Pilsen, Czech Republic, 91–96 [15] HERNÁNDEZ-BERMEJO, B., FAIRÉN, V., BRENIG, L. Algebraic recasting of nonlinear ODEs into univer- sal formats, J. Phys. A, Math. Gen., 1998, 31, 2415– 2430 [16] HERNÁNDEZ-BERMEJO, B., FAIRÉN, V. Lotka- Volterra representation of general nonlinear sys- tems, Math. Biosci., 1997, 140, 1–32 [17] HERNÁNDEZ-BERMEJO, B. Stability conditions and Lyapunov functions for quasi-polynomial sys- tems, Applied Mathematics Letters, 2002, 15, 25– 28 [18] BOYD, S., GHAOUI, L. E. , FERON, E., BALAKR- ISHNAN, V. Linear Matrix Inequalities in System and Control Theory (SIAM, Philadelphia) 1994 [19] SCHERER, C., WEILAND, C. Linear Matrix Inequalities in Control (DISC, http://www.er.ele.tue.nl/sweiland/lmi.pdf), 2000 [20] KOCVARA, M., STINGL, M. A code for con- vex nonlinear and semidefinite programming, Opti- mization Methods and Software, 2003, 8, 317–333 [21] TUAN, H., APKARIAN, P., NAKASHIMA, Y. A new Lagrangian dual global optimization algorithm for solving bilinear matrix inequalities, International Journal of Robust and Nonlinear Control, 2000, 10, 561–578 [22] CAO, Y.-Y., LAM, J., SUN, Y.-X. Static Output Feedback Stabilization: An ILMI Approach, Auto- matica, 1998, 12, 1641–1645 [23] KOCVARA, M., STINGL, M. TOMLAB/PENBMI solver (Matlab Toolbox) 2005, PENOPT Gbr. [24] I. THE MATHWORKS: MATLAB environment 2007) copyright 1984-2007 The MathWorks, Inc. http://www.mathworks.com [25] HANGOS, K. M., SZEDERKÉNYI, G. The under- lying linear dynamics of some positive polynomial systems, Physics Letters A, 2012, 376(45), 3129– 3134 Advertise upcoming meetings, conferences, and workshops; make public announcements; introduce your research laboratory; a new product or a service in the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-298, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 59-64 (2013) A STRUCTURED MODEL-BASED DIAGNOSIS METHOD FOR DISCRETE DYNAMIC PROCESSES USING EVENT SEQUENCES ATTILA TÓTH B1, KATALIN M. HANGOS1,2, AND ÁGNES WERNER-STARK1 1Department of Electrical Engineering and Information Systems, University of Pannonia, Veszprém H-8200, HUNGARY 2Computer and Automation Research Institute, Budapest, HUNGARY B E-mail: atezs82@gmail.com A novel model-based fault detection and diagnosis method is proposed that is based on following event sequences mea- sured in a discrete dynamic process. The model of the nominal and faulty operation modes is given in the form of event sequences, that are decomposed according to the components and sub-components present in the process system. The faulty event sequences are defined using extended procedure HAZID tables. A diagnostic algorithm is also presented that uses a component-wise decomposed form of the event sequences. The operation of the algorithm is illustrated on a simple example of a process system consisting of three similar tanks. Keywords: Process monitoring, diagnosis, discrete event systems, qualitative models Introduction Fault prevention and mitigation in the field of process system management is a task of crucial importance in avoiding serious accidents. Thus, numerous hazard iden- tification (HAZID) techniques have been developed in the past decades to ensure the safe operation of process systems and to relieve effects of faults (see [1] for a broad presentation of the field). Among these techniques, the most important methodologies involve the function- driven HAZOP (HAZard and OPerability, see [2]) anal- ysis and the component-driven FMEA (Fault Mode and Effects Analysis). There have been results in the past decade for automating the creation of HAZOP analysis ([3] with a concrete application described in [4]). Blend- ing the component-driven and function-driven analyses also resulted in a novel hazard identification approach de- scribed in [5]. Although the information collected in the HAZOP and FMEA studies serve the purpose of hazard identifi- cation, these studies can be the basis of diagnostic pro- cedures, too. A model-based diagnostic method based on HAZOP and FMEA information is reported in [6]. It is important to note that the above techniques con- centrate on the static case when the deviation from a nor- mal steady-state behaviour is of importance. Therefore, the transient case when the plant is controlled by an op- erational procedure is not addressed in these results. A recent study [7] tries to deal with the diagnostic task by using a specially constructed P-HAZID analysis and a diagnostic algorithm. In this paper, this diagnostic idea is extended to be able to handle more complex diagnos- tic tasks - by taking advantage of a possible decomposi- tion of typical process systems along their similar com- ponents. Basic notions In a complex system the full dynamic model that de- scribes its behaviour under normal and faulty operation models is rarely available, therefore one should base the diagnosis on qualitative information both in terms of the dynamic models and in the measured data. Here we briefly summarize the basic notions for qualitative model based diagnosis. Qualitative range spaces Current values of continuous measurable outputs in pro- cess systems can be described using a properly selected qualitative range space. For example, to describe the value of a level sensor in a tank, the following range space can be used: Qe = {e�, 0, L,N,H, e+}. (1) Here, 0 means an empty tank, L, N and H means low, normal and high fluid level, respectively, while e� and e+ refer to unmeasurably low and high fluid levels (this might mean a failure in the level sensor itself). This range space will be used to describe system outputs during op- eration. Input-output event sequences Operational procedures in process systems are detailed list of instructions for the plant operator personnel to 60 perform certain operations on the plant. Procedures can be formally described using finite input-output event se- quences where a single event describes a change in ei- ther the inputs or the outputs of the system at a specific time instant. Therefore the syntax of a single input-output event (at time instant t) is the following eventt = (t; input values; output values). The input in an event always refers to a state of an actua- tor component in the process system (e.g., in the case of a valve it can be open (op) or closed (cl)). On the other hand, the output in an event refers to a value of an output of the process system in the qualitative range space using the qualitative set defined in Eq. (1). Sequences formed from these events are called traces and defined as T (t1, tn) = eventt1 , ..., eventtn . Separate events in a trace contain the same inputs and outputs. Note that the discrete event time instances ti are abbreviated by their indexes i in the description, i.e. 2 stands for t2. Examples of single events for a two-input single out- put case include (1; cl, op;N) , (2; cl, op; L) where at time instances t1 and t2 the two valve inputs are held closed and open, respectively, while the level de- creases form its normal value to the low level. The trace formed form the above two consecutive events is written as T (1, 2) = ((1; cl, op;N), (2; cl, op; L)) . For every operational procedure there exists a trace (called the nominal trace) which describes its behaviour under fault-free conditions. The diagnostic method com- pares this trace to other traces which may have been executed under faulty conditions (called characteristic traces), and the differences (called deviations) are later used to find possible malfunctions of components in the system. Deviations Nominal and characteristic traces can be compared by comparing their corresponding event fragments. The dif- ference between two corresponding event fragments is described by a deviation. Deviations are formed from a deviation guideword and the nominal event from which the corresponding characteristic trace event is deviating from. The following deviation types are used during di- agnosis: • never-happened: When the particular event never happened in the characteristic trace. • later: When the event happened in the characteristic trace, but at a later time instant. • earlier: When the event happened in the character- istic trace, but at an earlier time instant. • greater: When a particular output’s qualitative value was higher in the characteristic event than that of the nominal trace. • smaller: When a particular output’s qualitative value was lower in the characteristic event than that of the nominal trace. For the detailed description of the greater and smaller qualitative relations, please refer to [7]. Procedure HAZID As a combination and extension of the widely used FMEA and HAZOP analyses (for details, refer to [7] or [1], and in particular to [5]), the procedure HAZID (ab- breviated as P-HAZID) analysis can be used for fault di- agnosis during operational procedures in a given process system. The result of this P-HAZID analysis is given in the form of a spreadsheet and it consists of deviations to- gether with their implications and possible (root) causes. A cause is considered to be a root cause if it is a non- measurable failure mode of a system component (which is an elementary part of the system). For example, a leak on a tank is considered as a root cause. A simple example of a P-HAZID table can be found in Table 1. Using the initial set of differences (deviations) be- tween the characteristic trace and the nominal trace, the set of possible root causes can be found using simple rea- soning. For details, refer to Ref. [7]. The diagnostic algorithm uses this technique first to find possible P-HAZID row(s) to start from (using the set of initial deviations). Then, following the deviation chains defined by these rows, the algorithm proceeds to- wards a possible root cause by traversing new rows based on the initial set of deviations. Using this procedure, it may end up at a root cause or at a row with deviations from which it cannot proceed forward, because they are not contained in the initial set of deviations. The algo- rithm assumes that the root causes are static and they hap- pened before the execution of the procedure began. Component based diagnosis It is widely known that complex systems can often be decomposed in a hierarchical way using simple non- dividable elements that are called components. The con- nection of these components is usually specified in terms of a graph called flowsheet. Such a decomposition can be used effectively for the operation of a reasoning-based diagnostic algorithm. Component based structural decomposition The above mentioned fault diagnosis based on the P- HAZID analysis is only developed for process systems 61 Table 1: A simple example of a P-HAZID table. Inputs: op=open, cl=closed. Outputs: 0=no, L=low, N=normal. Deviations: NH=never-happened, LAT=later, EAR=eariler, SML=smaller and GRE=greater. Faults: TANK-LEAK is the leak of the tank and POS-BIAS is the positive bias failure of the tank level sensor. Cause Deviation Implication TANK-LEAK NH(2;op,cl;L) NH(3;op,cl;N) NH(2;op,cl;L) NH(3;op,cl;N) NH(4;op,op;N) TANK-LEAK SML(2;op,cl;L) SML(3;op,cl;N) SML(2;op,cl;L) SML(3;op,cl;N) SML(4;op,op;N) POS-BIAS GRE(1;op,cl;0) GRE(2;op,cl;L) GRE(1;op,cl;0) GRE(2;op,cl;L) GRE(3;op,cl;N) GRE(2;op,cl;L) GRE(3;op,cl;N) NH(4;op,op;N) POS-BIAS NH(1;op,cl;0) EAR(2;op,cl;L) NH(1;op,cl;0) EAR(2;op,cl;L) EAR(3;op,cl;N) consisting of different individual components in [7]; the possible redundancy of such systems (e.g. multiple com- ponents of the same characteristics) were not taken into account. However, complex process systems in practice can be decomposed into a connected network of more simple components. For example, the process system in Fig. 1 can be decomposed into three smaller similar com- ponents each formed by an input and an output valve and a tank. When developing the decomposed form of a process system it can happen that some elements are part of mul- tiple subsystems as in the case of valves VB and VC in Fig. 1. These elements are called boundary components, and are assumed to be error-free during the diagnosis. Traces affecting different components can also be de- composed into a chain of trace fragments each referring to a single component of the trace. Events in such a trace fragment have only a subset of inputs and outputs of the united trace (only the inputs and outputs of the particular component that is present in them). Fragments also have information about the next trace fragment (called the next trace), and there is a starting condition (an event) that is associated with them to help the diagnosis. Along with the trace fragment, each component has its own associ- ated P-HAZID spreadsheet. The component based diagnostic algorithm Applying the diagnostic approach described in [7] on a decomposed process system, the components can be diagnosed separately against faults, treating them as a whole system during diagnosis. After the separate diag- nosis, the root causes can be collected and the result- ing set of root causes yields to the set of root causes in the united system. Using the component decomposition, the size of the HAZID information required can be made lower in cases when similar connected subsystems form the process system to be diagnosed. On the other hand, VA VB VC VD TC TB TA TA LEVEL TB LEVEL TC LEVEL Figure 1: Process system consisting of 3 similar tanks system-level deviations need to be converted into compo- nent deviations by aligning their time and reducing the in- puts and outputs in their associated events to component- level inputs and outputs. The diagnosis is then performed by comparing the whole nominal trace with the characteristic trace, and then distributing the deviations (differences) among the components. Before distributing, time alignment and re- duction of input and output states to component level are performed (including the component’s boundary el- ements). After the distribution, component-level diagno- sis may begin to explore possible (root) causes on the component-level. The diagnosis proceeds along the flow direction by starting from the first component, takes the deviations, generates the set of possible root causes from them, and then tries to proceed to the next component by checking the starting condition of the next trace fragment - if there is any. If the start condition is fulfilled, the diagnosis con- tinues, otherwise it halts. For example, in the case of Fig. 1 the second fragment might have a start condition con- taining a statement about the minimum level of fluid in tank TA, and in the case of the congestion of valve VA no fluid is coming into a system, therefore, tank TA is not even filling up to the specified minimum level. In this case the diagnosis stops. The result of the diagnostic algorithm is always the union of identified and non-identified root causes created by the consecutive diagnostic algorithm that runs on the components of the consequent nominal trace fragments. For reference, the whole diagnostic algorithm is pre- sented as a pseudo-code in Algorithm 1. The algorithm collects all root causes (sets INC and IRC) given a com- ponent decomposition, a starting component, and a pos- sibly faulty characteristic trace. Case study In the case study the diagnosis procedure for the sim- ple process system in Fig. 1 containing three sequentially connected tanks is used. The measured output of the tanks 62 Algorithm 1 Component-based reasoning procedure 1: INC {;} 2: IRC {;} 3: actualComponent startComponent 4: continue true 5: shift 0 6: while continue do 7: DEV GENERATEDEVIATIONS(actualComponent, chrTrace, shift) 8: FDP COLLECTFINALDEVIATIONPAIRS(DEV) 9: for all pair 2 FDP do 10: startDeviation proj1(pair) 11: startImplication proj2(pair) 12: REASON(startDeviation, startImplication, actualComponent.phazid) 13: end for 14: if actualComponent has successive component and its start condition evaluates to true then 15: shift length(actualComponent.trace)� 1 + shift 16: actualComponent GETCOMPONENT(actualComponent.successiveComponent) 17: else 18: continue false 19: end if 20: end while 21: function GENERATEDEVIATIONS(component, chrTrace, shift) 22: DEV {;} 23: nomTrace component.trace 24: reducedTrace TRIMTRACE(chrTrace, shift, shift+ component.trace.length) 25: locChrTrace CONVERTTOCOMPONENTLEVEL(reducedTrace) 26: for T := 1 to length(nomTrace) do 27: for all deviation D of locChrTrace from nomTrace at time T do 28: DEV DEV S (D) 29: end for 30: end for 31: return DEV 32: end function 33: procedure REASON(deviation, implication, phazid) 34: if 9R 2 ROWS(phazid), deviation = devphazid(R), implication = impphazid(R) then 35: for all {R, devphazid(R) = deviation and impphazid(R) = implication} do 36: if causephazid(R) 2 RC then 37: IRC IRC S causephazid(R) 38: return 39: else if causephazid(R) 2 DEV and causephazid(R) � devphazid(R) in DEV then 40: REASON(causephazid(R), devphazid(R), phazid) 41: else 42: INC INC S causephazid(R) 43: return 44: end if 45: end for 46: else 47: INC INC S causephazid(R) 48: return 49: end if 50: end procedure 63 Table 2: Tank fill operational procedure. Time Input values Output values VA VB VC VD TA TB TC 1 op cl cl cl 0 0 0 2 op cl cl cl L 0 0 3 op op cl cl N 0 0 4 op op cl cl N L 0 5 op op op cl N N 0 6 op op op cl N N L 7 op op op op N N N Table 3: Normal fill in a single tank with no faults. Input valve Output Valve Tank Level op cl 0 op cl L op op N is the tank level that takes its values from Qe in Eq. (1), and the input variables are the valve positions (open (op) or closed (cl)). Components, operational procedure, and nominal trace Every tank may contain no fluid (the tank level is equal to no (0)), may be low on fluid (level value is low (L)), or might have normal fluid level (level value is normal (N)). In every time instant the level increases by one qualitative magnitude (i.e. from no to low or from low to normal) if fluid is coming through the input valve but the output valve is closed. Due to the same size of the valves the effect of fluid flow out of the system is similar, but in the opposite direction (from normal to low or from low to no). The valve positions (open (op) or closed (cl)) can be changed by the operator; they are considered as inputs of the system. Leak in the tank is assumed to be equal to the size of an open valve (i.e. a quite substantial leak). The considered operational procedure is the initial filling of all the three tanks with fluid (the fill operational procedure in short), and is described in detail in Table 2. The process system can be decomposed into three components, therefore the fill operational procedure can also be partitioned into three identical procedure frag- ments along the component boundaries (the VB and VC valves). The fragment abstracted from the three identical trace fragments associated to the three tanks can be observed in Table 3. It has only the subset of inputs and outputs which are directly related to the particular tank component - the input and output valve and the tank level. The corresponding component P-HAZID table can be found in Table 4 with some of the component faults and deviations associated to them. Instances of this P-HAZID table are used in the case of all three tanks during diag- nosis. Table 4: P-HAZID table of a single tank component with two valves for a reference trace of Table 3. Faults: TANK-LEAK is leak of the tank, POS-BIAS is the pos- itive bias fault of the level sensor and NEG-BIAS is the negative bias failure of the level sensor. Cause Deviation Implication TANK-LEAK NH(2;op,cl;L) NH(3;op,op;N) TANK-LEAK SML(2;op,cl,L) SML(3;op,op;N) NEG-BIAS LAT(1;op,cl;0) NH(2;op,cl;L) LAT(1;op,cl;0) NH(2;op,cl;L) NH(3;op,op;N) NEG-BIAS SML(1;op,cl;0) SML(2;op,cl;L) SML(1;op,cl;0) SML(2;op,cl;L) SML(3;op,op;N) POS-BIAS NH(1;op,cl;0) EAR(2;op,cl;L) NH(1;op,cl;0) EAR(2;op,cl;L) NH(3;op,op;N) POS-BIAS GRE(1;op,cl;0) GRE(2;op,cl;L) GRE(1;op,cl;0) GRE(2;op,cl;L) GRE(3;op,op;N) Table 5: Tank fill operational procedure with a leak in the second tank TB. The leak caused two different events in the operational procedure related to TB (in bold), these differences resulted in the four deviations the diagnosis could start from. Time Input values Output values VA VB VC VD TA TB TC 1 op cl cl cl 0 0 0 2 op cl cl cl L 0 0 3 op op cl cl N 0 0 4 op op cl cl N 0 0 5 op op op cl N 0 0 6 op op op cl N 0 0 7 op op op op N 0 0 The operation of the diagnostic algorithm The operation of the diagnostic algorithm is illustrated with the case when a rupture of the second tank is present as a root case (fault). A characteristic trace with the leak of the second tank can be seen in Table 5. The size of the leak is assumed to be larger or equal to the size of an outbound pipe, there- fore, the tank cannot fill up and no fluid can flow to the third tank TC. The starting condition of the TB and TC tank compo- nents is the appropriate “normal” level in the preceding tank. In that way it is ensured that diagnosis is done on the operational components only. The diagnosis of this faulty scenario begins by start- ing with the first tank component TA. There are no differ- ences (and therefore no deviations) regarding this com- ponent. The start condition of the second component is fulfilled, therefore, the diagnosis moves towards the next component TB. However, in the case of TB the following deviations are found after comparing the nominal and characteristic traces (due to that two consequent events did not hap- 64 pened, instead, two events happened with lower output values at time instant 4 and 5 in the operational proce- dure): • never-happened(4;open,closed;low) • never-happened(5;open,open;normal) • smaller(4;open,closed;low) • smaller(5;open,open;normal) The time instant of the deviations are shifted back by 2 units because the second component’s first event happens at the third system-level time instant. After that, the diag- nosis is initiated on the HAZID table. Searching for the already found deviations and connecting them to possible root causes (as in the case of the original diagnostic idea in [7]) the leak of the second tank can be found. Because of the lack of fluid in the second tank, the start condition of the third component is not fulfilled, therefore, the diagnostic process halts at this step result- ing in a single possible root cause being the TANK- LEAK. Conclusions A novel component based extension of the single compo- nent diagnostic algorithm presented in [7] is described in this paper. Using the extension, the domain of application can be extended to more complex composite process sys- tems. Driven by the decomposition of the overall system into components, the P-HAZID tables used for diagnosis are processed at component level by the diagnostic algo- rithm. The extended method is efficient in the cases when the overall process system consists of similar small com- ponents. The component-based diagnostic procedure was de- scribed on a formal level, along with its proposed pseudo- code. A case study for a process system of multiple com- ponents and a simple failure was also provided. The following improvements are planned to extend the component-based diagnostic approach: • The procedure is based on the assumption that the boundary elements between different components are free of failures. As a future work, this limitation might be removed by using a higher level reasoning above the components (as in the form of a system- level HAZID table, for example). • At the moment, the algorithm is only working for al- ready coded static event information in order to find faults in the system. Diagnosis would be more valu- able if events could be processed dynamically, thus the diagnostic procedure could be executed real- time along with the operational procedures. • Diagnosis would be more accurate if the derivatives of internal states (e.g. the derivative of the tank level) were present in the events. Acknowledgements The research was supported by the Hungarian Research Fund through grant K83440. We acknowledge the fi- nancial support of the Hungarian State and the Eu- ropean Union under the TAMOP-4.2.2.A-11/1/KONV- 2012-0072. REFERENCES [1] CAMERON, I. T., RAMAN, R. Process Systems Risk Management. Vol. 6 of Process Systems Engineering (Elsevier Academic Press, San Diego, CA) 2005 [2] AS IEC 61882-2003: Hazard and operability studies (HAZOP studies) Application Guide [3] VENKATASUBRAMANIAN, V., RENGASWAMY, R., YIN, K., KAVURI, S. N. A review of process fault detection and diagnosis Part I: Qualitative model- based methods, Computers and Chemical Engineer- ing, 2003, 27(3), 293–311 [4] VENKATASUBRAMANIAN, V., ZHAO, J. S., VISWANATHAN, S. Intelligent systems for HAZOP analysis of complex process plants, Computers and Chemical Engineering, 2000, 24(9-10), 2291–2302 [5] SELIGMANN, B. J., NÉMETH, E., HANGOS, K. M., CAMERON, I. T. A blended hazard identification methodology to support process diagnosis, Journal of Loss Prevention in the Process Industries, 2012, 25 , 746–759 [6] NÉMETH, E., LAKNER, R., CAMERON, I., HAN- GOS, K. M. Fault diagnosis based on hazard iden- tification results, in Preprints of the 7th IFAC Sym- posium on Fault Detection, Supervision and Safety of Technical Processes, Barcelona, Spain, June 30 - July 3, pp. 1515–1520 [7] TÓTH, A., HANGOS, K. M. , WERNER-STARK, A. HAZID information based operational procedure di- agnosis method, in 12th International PhD Workshop on Systems and Control, Veszprém, 2012. aug. 27 (A. MAGYAR, ed.), ISBN 978-615-5044-71-7 (Uni- versity of Pannonia, Veszprem), pp. 1–6 (on CD) HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 65-75 (2013) ROLE OF STEADY STATE DATA RECONCILIATION IN PROCESS MODEL DEVELOPMENT BARBARA FARSANGB1, SÁNDOR NÉMETH1, AND JÁNOS ABONYI1 1University of Pannonia, Department of Process Engineering, H-8200 Veszprém, Egyetem Street 10., HUNGARY BE-mail: farsangb@fmt.uni-pannon.hu In chemical and hydrocarbon industry operational efficiency is improved by model-based solutions. Historical process data plays an important role in the identification and verification of models utilized by these tools. Since most of the used information are measured values, they are affected by errors influencing the quality of these models. Data reconciliaton aims the reduction of random errors to enhance the quality of data used for model development resulting in more reliable process simulators. This concept is applied to the development and validation of the complex process model and simulator of an industrial hydrogenation system. The results show the applicability of the proposed scheme in industrial environment. Keywords: data reconciliation; flowsheeting simulator; validation; balance equation; industrial hydrogenation system Introduction In recent years, development of industrial technologies has been determined by the proliferation of computing and information technology. Today’s technologies are characterized by widespread application of process engineering tools. The success of these tasks depends on the accuracy of data and applied models [1]. Thanks to the evolution of information technology, on-line and historical process data — coming from chemical process systems — are available. The collected data provide the opportunity for engineers to better understand the processes, anomalies, and malfunctions [2]. Monitoring of process variables allows us to ensure the consistent product quality. Collected data can also be used for the development and validation of process simulators. However, measurements are always affected by errors during the measurement, processing and transmission of the measured signal. Errors (gross, random, bias) in measured data affect the quality of process models and can lead to significant decrease in plant performance. Estimation of true conditions of process states is important to achieve optimal process monitoring, control and optimization. Therefore, several methods have been developed to minimize measurement errors thereby enhancing the reliability, accuracy and precision of data. First, analogue and digital filters were used to reduce the effect of high frequency noise [3]. Gross errors were detected with various data validation methods, which include checking whether the measured data and the rate at which it is changing is within predefined operational limits. Nowadays, smart sensors are used for determine whether there is any hardware problem or the measured data are appropriate. New methods were statistical quality control tests that are applied to each measured variable separately. Although these methods improve the reliability and accuracy of the measured data, they do not ensure consistency of the data with respect to the inter-relationships between different process variables. Therefore new methods – data reconciliation and gross error detection – are developed in chemical engineering [4]. The main difference between data reconciliation and other filtering methods is that data reconciliation uses process model constraints and that the results satisfy constraints and balance equations. Data reconciliation techniques can be applied to reduce random errors of measurements, while the other techniques mentioned above do not. Gross error detection can be used for eliminate systematic errors so simultaneous data reconciliation and gross error detection have emerged as a key of online optimization [5] (see Fig. 1). Data reconciliation techniques take minimal corrections of the measured variables to satisfy a set of model constraints. Based on the difference between the measured and reconciled data, the following questions can be answered. Can we consider a related set of measurements acceptable based on our previous knowledge of the system? Are the measurements consistent? If not, what can be the source of the error? Based on the available measurements and a prior knowledge, what is the most likely state of the system? Answers to these questions are important in development of technology (e.g., monitoring, optimization, simulation, control, instrument maintenance). Using this technique, we can verify the acceptability of measurements, improve the accuracy of measurements, estimate model parameters and unmeasured variables, and it can be used 66 Applications S ys te m s ta tu s id en tif ic at io n R e a l ti m e d a ta b a se Reconciliation Gross error detection Variable classification Steady state identification Data filtration Detection of abnormalities Figure 1: General methodology for on-line data correction [6] for fault detection. In order to ensure the consistency of measured process variables, models are used in data reconciliation. During data reconciliation, information is obtained from both measurements and process models. Depending on the types of the models, data reconciliation techniques can be separated into two different problems: steady state and dynamic (see Fig. 2). In the case of steady state data reconciliation, model constraints are algebraic equations. If we are dealing with dynamic processes, differential equations are used as model constraints. In the case of linear data reconciliation, model constraints are defined by linear equations. The simplest data reconciliation technique is the steady state linear method that can describe, e.g., the mass balance of technology. Moreover, most chemical processes have nonlinear characteristics and constraints are nonlinear equations; this problem can be called as nonlinear data reconciliation. The goal of the presented research is to propose a methodology to support the development and validation of complex process models and simulators based on increasing the quality of measurement data used for (kinetic) parameter identification and validation. In section Literature survey, a historical background of steady data data reconciliation is briefly reviewed. The basic method of linear and nonlinear steady state data reconciliation is described in section Theoretical background of steady state data reconciliation. Data reconciliation cannot be used in process development only, but it is a useful tool in model improvement too. The principle of this application is introduced in section Combined application of data reconciliation and flowsheeting simulator. Our approach aims the reduction of random errors to enhance the quality of data using data reconciliation and flowsheeting simulator simultaneously. The method is described in subsection The developed model based data reconciliation technique. The proposed approach is illustrated based on an industrial hydrogenation system. Figure 2: Types of data reconciliation [7] In section Case study, the analyzed technology is introduced and the applicability of the method is illustrated on the basis of the mass balance of the hydrogenation system. Then, some consequences are drawn in section Conclusion. Literature survey KUHAN and DAVIDSON presented the first study about minimizing measurement errors in chemical process industry using a data reconciliation technique in 1961 [4]. They described the steady state linear data reconciliation technique and demonstrated a general solution for mass balance of technology. This seminal paper has started a new research area in process engineering. NOGITA [8] and MAH [9] improved a new method capable of detection of gross errors. Later, MAH [10] found that the base method requires a process system in steady state conditions and that the process constraints must be linear. CAO and RHINEHART confirmed that the methodology will be successful if process data come from steady state conditions [11]. In these cases, data reconciliation problems can be solved with standard numerical methods like sequential quadratic procedures [12] and quasi-Newton algorithms [13]. A multiple Gauss-Jordan elimination algorithm was published by MADRON and VEVERKA [14]. The classification method is based on matrix decomposition and permutation. The linear method was improved by CROWE (1986) who published a new technique to solve nonlinear data reconciliation using matrix projection [15]. Furthermore, BAGAJEWICZ and JIANG [16] presented an improved method that is applicable for dynamics system. The developed method from PRATA [17] is suitable for examination of dynamic nonlinear process systems. Despite of the numerous good reviews of data reconciliation techniques, industrial application of this method started only in the 1980s. There are several complex software packages supporting data 67 Table 1: Industrial application of data reconciliation (newer industrial scenario and publications are added to PRATA’s summary [18]) Industrial scenario Author (year) Absorption refrigeration systems MARADIAGA et al. (2013) [19] Beverage alcohol distillation plant MEYER et al. (1993), NOUNOU and BAKSHI (1999), SCHLADT and HU (2007) Chemical extraction plant HOLLY et al. (1989) Ethylene and ammonia plant SANCHEZ et al. (1992), SANCHEZ and ROMAGNOLI (1996), PLÁCIDO and LOUREIRO (1998) Exxon chemical process MCBRAYER et al. (1998) , SODERSTROM et al. (2000) Gas pipeline systems BAGAJEWICZ and CABRERA (2003) Gases network in an iron and steel making plant YI and HAN (2004) Hydrogen plant BUSSANI et al. (1995), CHIARI (1997), SARABIA et al. (2012) [20] Industrial coke-oven-gas purification process FABER et al. (2006) HU and SHAO (2006) Industrial distillation column ISLAM (1994), WEISS et al. (1996), SANCHEZ et al. (1996), BOUROUIS et al. (1998), LI et al. (2001), BHAT and SARAF (2004), CHATTERJEE and SARAF (2004), CHEN et al. (2013) [21], KELLER et al. (2012) [22] Industrial ETBE reactor DOMINGUES et al. (2012) [23] Industrial furnace PIERUCCI et al. (1996), EKSTEEN et al. (2002) Industrial hydrometallurgical plants for a gold extraction DE ANDRADELIMA (2006) Industrial polymerization reactor VIEIRA et al. (2003), PRATA et al. (2006, 2008, 2009, 2010 [24]) Industrial synthesis gas for production of ammonia CHRISTIANSEN et al. (1997) Industrial utility plant LEE et al. (1998) Methyl-terc-butyl-ether plant AL-ARFAJ (2006) Mineral and metallurgical plants VASEBI et al(2012) [25] Nuclear power reactor VALDETARO (2011) [26] Refinery PICCOLO and DOUGLAS (1996), ZHANG et al. (2001) Sulfuric acid plant CHEN et al. (1998), OZYURT and PIKE (2004) Turbine cycle of a boiling water reactor SUNDE and BERG (2003) Vinyl acetate and ketene plants DEMPF and LIST (1998) Water processes MARTINS et al. (2010) [27] 68 reconciliation (for instance SigmaïňĄne, DataCon, Vali, Inlibra). These tools have interfaces to information systems used to manage process data and have special model building functionalities that allow almost automatic building of plant-wide balance equations. Nowadays, data reconciliation is widely applied in various processing industries (e.g. refining, chemical industry, metals, mining, and power industry). Despite of this fact, relatively few articles deal with industrial applications (PRATA et al. prepared a summary in 2009 [18]; newer applications and papers are collected in Table 1). STANLEY et al. classified variables as observable and unobservable [28]. CROWE categorized the variables from another point of view: the basis of classification is the given measured variable can be calculated from other measured variables using process models or not. The names of the two classes are redundant and non- redundant measurements [29]. ALI and NARASIMHAN applied graph theory to analyse the sensor network, classify the variables, and compute the redundancy degree of each variable. They claim that not only measured variables can be redundant. Those unmeasured variables that can be estimated in multiple ways are also termed as redundant [30]. Theoretical background of steady state data reconciliation In this section, the basis of linear and nonlinear data reconciliation techniques are described. Different values are distinguished. Real values (without error) of variables are designated by diacritic caron (̌ ), measured values by tilde (̃ ), while estimated values by circumflex (̂ ). In general, optimal estimates for process variables by data reconciliation are solutions to a constrained least-squares or maximum likelihood objective function, where measurement errors are minimized with process model constraints. The steady state data reconciliation problem can be formulated as an optimization problem by minimizing J(y̌, ž) = (ỹ � y̌)T V�1 d̃ (ỹ � y̌) (1) subject to f(y̌, ž) = 0 (2) and g(y̌, ž) � 0. (3) The solution of the optimization problem is performed with the following simplifying assumptions • The measurement error is independent from the balance variables. • The expected value of measurement error is zero. • The measurement error is in a normal (Gauss) distribution. • Errors of measurements are independent from each other (diagonal covariance matrix). • The covariance matrix is positive definite, so its inverse matrix exists. In the following two subsections, solution methods are presented in the cases of linear and nonlinear models. Linear data reconciliation In the linear case, when all variables are measured, Eq. (1) means minimizing J(y̌, ž) = (ỹ � y̌)T V�1 d̃ (ỹ � y̌) subject to Ay̌ = 0. (4) Gauss distribution (d̃ ⇠ N � 0,Vd̃ � ) is assumed, so the density function of measurement error is f ⇣ d̃ ⌘ = C · exp h �0.5 · d̃TV�1 d̃ d̃ i (5) from which it follows that f (ỹ) = C · exp h �0.5 · (ỹ � y̌)T V�1 d̃ (ỹ � y̌) i (6) and that the f̌ = Ay̌ � b = 0 (7) criterion is satisfied in accordance with the balance equation. The essence of the most likely estimation: if ỹ is given and its density function is known, which ỹ ! ŷ parameter of density function will be f (function vector of equality model constraints). Since the logarithm function is strictly monotonously increasing, the maximum of the logarithm of the f(y) function recording only positive values is at the same place where the maximum of the argument is. Exploiting it on the above function and using the Lagrange multipliers method, the result is a vector equation system with two unknowns: V�1 d̃ ŷ �AT �̂ = V�1 d̃ ỹ (8) Aŷ = b (9) The solution of the vector equation system is ŷ = ⇣ I�Vd̃A T � AVd̃A T ��1 A ⌘ ỹ + Vd̃A T � AVd̃A T ��1 b (10) Thus, the estimation requires a matrix-vector multiplication and a vector addition. Because ỹ is only the variable, in the other operations there are only constants [31]. In practice, not all streams and properties are measured in plant due to physical and economical reasons. The data reconciliation technique is suitable for estimating unmeasured variables. In this case, the optimization problem can be solved by the method of 69 Projection Matrix. In the first step, the incidence matrix is separated into two parts: Ayy̌ + Az ž = 0. (11) Unmeasured flows (ž) should be eliminated by pre- multiplying both sides by a projection matrix, such as PAz ž = 0. Eq. (4) then can be rewritten as min J(y̌, ž) = (ỹ � y̌)T V�1 d̃ (ỹ � y̌) (12) subject to PAyy̌ = 0. (13) The projection matrix is from Q-R factorization of matrix Az: Az = QR = [Q1 Q2] [R1 0] T (14) The projection matrix is the transformation of Q2 matrix (P = QT 2 ). If A matrix is replaced by matrix PAy in Eq. (10), then ŷ = ⇣ I�Vd̃(PA)T � (PA)Vd̃(PA)T ��1 (PA) ⌘ ỹ + Vd̃(PA)T � (PA)Vd̃(PA)T ��1 b. (15) After reconciled values of measured variables are obtained, the next step is to estimate the unmeasured variables using the information from process models (Eq. (11)) as Ayy̌ = �Az ž (16) This linear equation can be solved to get estimated values of the unmeasured variables. Usually, the number of equations is greater than the number of unmeasured variables. The least-squares technique can be applied and the solution is ẑ = �(AT z Az) �1AT z (Ayŷ). (17) Nonlinear data reconciliation The data reconciliation problem can be extended to nonlinear steady state models. It is necessary because processes in chemical industry cannot be described only by linear models; behaviour of chemical processes are often nonlinear. If we wish to simultaneously reconcile mass flow, composition, or temperature measurements, mass, component, or energy balances have to be included as constraints. Sometimes, inequality constraints have to be defined (e.g. mass/mole fraction of components have to be in the [0:1] interval). Moreover, if we take thermodynamic equilibrium relationships and complex correlations for thermodynamic and physical properties as constraints, nonlinear data reconciliation techniques have to be used. If all variables are measured, the bilinear data reconciliation problem can be reduced to a linear problem by introducing the “measured” compound flows (‘component flows’ or ‘energy flows’). Therefore, solution of the nonlinear data reconciliation problem is Figure 3: Integrated application of statistical and white box models [32] presented on the basis of the situation when measured and unmeasured variables occur. The objective function is the same as in the previous case (Eq. (1)). The difference is that constraints are nonlinear equations. This problem can be reduced to the linear case. This method is referred to as successive linearisation: nonlinear constraints can be linearised around working points using first-order Taylor series. The solution of the linearisation can be written as h = Byy̌ +Bz ž� f(y̌i, ži). (18) The initial values are usually the raw values. After linearisation the procedures of data reconciliation is described in subsection Linear data reconciliation (Q- R factorization, reconciliation of measured values, then estimation of unmeasured variables). If the reconciled and initial values are far from each other, new iteration begins. The initial value will be the result of the previous iteration. Different criteria can be defined: number of iteration or ||ŷn � ŷn�1|| is smaller than the specified tolerance values. Successive linearisation is a relatively simple and fast solution, but variable bounds cannot be handled with this method. Another option to solve nonlinear data reconciliation problem is using Nonlinear Programming (NLP) techniques that can estimate measured and unmeasured variables simultaneously. Sequential quadratic programming (SQP) and generalized reduced gradient (GRG) are usual techniques in handling nonlinear problems. These methods are more computationally demanding, but they are numerically more robust [33]. 70 Figure 4: Iterative model development Combined application of data reconciliation and flowsheeting simulator Mathematical models describe the connection between the relevant properties of the analyzed phenomena, process, or activities. In the case of complex chemical processes, mathematical models usually contain differential equations whose solution is not always possible analytically. The process simulator (it is often applied to improve complex industrial processes) contains models of the equipments so their mathematical description is unnecessary for the users. In addition, it includes different thermodynamic models and component databases. The flowsheeting simulator can be defined as “Use of a computer program to quantitatively model characteristic equations of a chemical process” [34]. The simulator is used in batch processing, integrated process engineering (for example, economic analysis and supply chain forecast, supervisory process control, on-line modeling and optimization, safety and reliability analysis), and process synthesis and design (for instance, heat integration, conceptual design). However, lots of information are needed from technology for the simulator to describe the real process exactly. Many of the needed data are measured (e.g., flow, composition, temperature, pressure). If the input of the simulator is faulty, the simulator cannot give reliable results. A new direction is the combination of simulators and data based statistical models that can handle random and gross errors [23, 35]. The structure of these hybrid models are shown in Fig. 3. • In the first case, the statistical model is the input of the physical model in the form of differential or algebraic equations or a complex flowsheeting simulator. In this case, the statistical model is used to estimate parameters and phenomena that are Figure 5: Options of data processing in flowsheeting simulator difficult to model. • Combination 2 shows the case when outputs of the physical model are transformed by a statistical model. • In the third option, the difference between measured and calculated variables are the inputs of the statistical model used for correction. Model-based data reconciliation techniques are similar to the third approach. Using data reconciliation technique random errors can be filtered. When accurate data are available, the flowsheeting simulator gives more reliable results that make further development of models possible. This means that data reconciliation allows us to check the reliability of these measurements [36] and reconsolidated data can be used to build accurate models. The developed model based data reconciliation technique Our goal is to design an expert system that can be used to check the acceptability of measurements and improve the flowsheeting simulator of technology. The proposed method is based on model-based data reconciliation techniques. Using data reconciliation, faulty measurements can be found. Random errors can be filtered using discrete Fourier transform. The frequency of the noise is much higher than the frequency of basic process so the discrete Fourier transform gives a graphical view about the differences between the basic process and the noise. Thus, the filtration of the frequency can be determined so the noise can be separated from the basic process. If the amplitude of the noise is known, the standard deviation can be estimated. Thereby, input of simulator is error-free, so systematic model mismatch can be recognized and the tuning of the model can be initiated. Data reconciliation requires an accurate model whose parameters require reconciled process data. In th case of an efficient iterative procedure, a properly working simulator can be developed (see Fig. 4). Since the simulator is able to process historical measured data, the difference between the theoretically achievable values and the measured outputs can give useful information (see Fig. 5). Four different output 71 Figure 6: Flowsheet diagram of the C3 hydrogenation system values can be compared in the case of the same variable. Measured output values are known. If all necessary information is available for data reconciliation, we get the reconciled output value. In addition, if process model is prepared in the flowsheeting simulator, we have two options for analysis: the input of the simulator can be the measured or the reconciled value. The difference between outputs can be used for two purposes. If the difference is insignificant, the measured data is acceptable and measuring instruments operate properly; calibration is not necessary. Moreover, it helps in the development of the flowsheeting simulator. If the calculated results with reconciled input are far from the reconciled output, the simulator do not describe the real process properly; maybe a parameter is not accurate or it needs structural change. The proposed method is illustrated on the basis of an industrial hydrogenation system. The case study shows examples for both cases; when the difference indicates the conformity of measurements and when it warns that some model parameters are incorrect. Case study The Tisza Chemical Group Plc. (TVK) is the largest petrochemical company of Hungary where polymer raw materials (ethylene, propylene, butylenes, etc.) are produced by steam cracking of naphtha or gasoline. At high temperature, numerous free radical reactions occur. Cracked gas includes light components that are produced in larger amounts, for example, methane, ethane, ethylene, acetylene, propane, propylene, methyl-acetylene, propadiene, n-butane, isobutane, 1,3-butadiene, etc. Some of them are undesirable, because they poison the catalysts for polymerization reactions (e.g. acetylene, methyl-acetylene, propadiene). Another problem is that separation of these components is difficult from the main products by distillation. Due to these problems, undesirable hydrocarbons are hydrogenated. Hydrogenation process of methyl- acetylene and propadiene is presented in this case study (see Fig. 6). The C3-selective hydrogenation process transforms methyl-acetylene and propadiene to propylene by catalytic reaction avoiding the transformation of propylene to propane. The methyl-acetylene and propadiene content of hydrogenated C3 fraction should be less than 1000 ppm. The concentration of the main components (methyl-acetylene, propane) is measured by online analysers (A1, A2 and A3) and there are two places where sample is taken twice a day for laboratory analysis (L1 and L3). Hydrogen and liquid olefin stream are fed to the reactor where chemical reactions take place. Outlet stream of the reactor is cooled by water so the gaseous C3 component condenses and the hydrogen and the liquid olefin phases are separated. The recirculated stream has two important functions: cools the reactor and dilutes the inlet C3 stream. Blowdown is required due to the accumulation of inert components. The most frequent type of hydrogenation reactors is the trickle-bed reactor in the olefin plant. Liquid olefin feed and gaseous hydrogen pass through the catalyst bed in the same direction: from top to bottom. Cooling of reactor is provided by vaporization liquid flow. The reactions occur in the top region of the reactor. The 72 Figure 7: Measured and filtered values of product C3 fraction reactions release a large amount of heat. The heat generation is influenced by the methyl-acetylene and propadiene (MAPD) content of the inlet C3 and by the partial pressure of the hydrogen. The MAPD content of the inlet flow of system is high, so the inlet flow of the reactor is diluted. In the range of 10 � 80 �C temperature, the average reaction heat of main reactions are the following: MA+H2 = C3H6 �Hr = �164.8 kJ/mole (19) PD+H2 = C3H6 �Hr = �172.2 kJ/mole (20) C3H6 +H2 = C3H8 �Hr = �124.6 kJ/mole (21) In the following sections, released heat vaporizes a part of the liquid phase. This phenomenon also ensures the cooling of the product flow. The simulator has been developed in Aspen Plus software as part of an expert system used for the monitoring and qualification of the operation of the technology. Calculations provide useful information about unmeasured variables, validity of on-line analyzers, and efficiency of catalysts. The residence time is small (less than one minute). The sample time was five minutes, so the steady state simulation gives the same result as the dynamic simulation. Historical process data entered the simulator by Aspen Simulation Workbook. The developed Excel and Visual Basic macro based framework allows the comparison of the measured data and the calculated results. Results and discussion The process has two inlet streams (C3 fraction and hydrogen) and two outlet streams (hydrogenated liquid C3 fraction and blowdown). Mass flows of every stream are measured. Based on the law of conservation of mass (there is no accumulation): Inlet C3 + H2 = Blowdown + Product C3 (22) Figure 8: Measured and reconciled values of product C3 fraction If we check whether the measured data satisfy the balance equation, we find that Eq. (22) is not fulfilled, so data reconciliation is necessary. Steps of the procedure are introduced based on mass flow of C3 product fraction. First time random error is separated from basic process using discrete Fourier transformation. Fig. 7 shows the measured and the filtered values of product C3 fraction. We defined the random noise as the difference between the measured and filtered values so the standard deviation and covariance matrix of random error (Vd̃) can be calculated. Since there is no accumulation (i.e., b = 0) Eq. (10) is simplified to ŷ = ⇣ I�Vd̃A T � AVd̃A T ��1 A ⌘ ỹ. (23) Since two input and two output streams are in the system, the incidence matrix is A = [1 1 � 1 � 1]. Every needed information are known so the matrix operation can be performed. Real industrial data are analyzed, so dimensionless units are included in the figures. The result of reconciliation is visible in Fig. 8. The result shows that the proposed method gives a minor improvement of the process values. Thanks to the validated accuracy of the mass measurements, these process values can be directly used for the validation of the simulator. In the next step, we compare the reconciled output with the values calculated by the simulator in two different ways. First, the input of the simulator is the measured data, then, the reconciled input data. The result of the comparison is shown in Fig. 9. First, we thought that the three curves coincide. However, as a portion is zoomed, the difference became apparent (see Fig. 10). If the input of simulator is the reconciled data, the calculated and the reconciled output are nearly the same. Although the difference between the two calculated curves may seem small in dimensionless space, do not forget that the small percentage difference can mean significant mass flow values in a real industrial process. 73 Figure 9: Calculated and reconciled data of product C3 fraction if the input of the simulator is the measured and reconciled data Figure 10: Reconciled and calculated outputs of simulator in case of product C3 fraction We calculated the square error in both cases (deviation from the reconciled output) and found that the difference decreased by 87% if the reconciled values are the input of the simulator. This shows that the simulator gives more accurate results if previously data are reconciled. The simulator supported by data reconciliation can be applied in fault diagnosis. Fig. 11 shows how the hydrogen flow changes over time. The continuous line represents the simulator results without data reconciliation, while the dashed line shows the reconciled values. Fig. 11 shows the similarity of the calculated and the measured values. After 300 hours of operation the calculated values exceed the measured data. There are several possible reasons. Less hydrogen goes to this system if • The MAPD content of the inlet stream is lower, because the mass flow of hydrogen is controlled in proportion of the MAPD mole flow of the input C3 Figure 11: Calculated and reconciled dimensionless mass flow of H2 if the input of the simulator is the measured data fraction. • After the regeneration when catalyst is active. • Other reactor is regenerated or reactivated. In this case, another reactor in the Olefin 2 plant was reactivated (after 300 hours). The reactivation process needed a lot of hydrogen so the hydrogen for the studied reactor was reduced. The received data did not include the rate of reduction. The result of the simulation showed that something was changed in the technology. This shows the technique can be used for monitoring and diagnostics of complex processes because it can be discovered if the operating condition has changed. Conclusion Data reconciliation is an effective method to obtain accurate process data that satisfy balance equations using process constraints. The application of steady-state data reconciliation is widespread in the process industry. Simulators are often applied to improve industrial processes, optimize operation, and identify bottlenecks of technology. Historical process data can be used for the identification and verification of models utilized by these tools. Usually, measured data do not satisfy balance equations, because all measurements are incorrect to some extent. Thus, it is necessary to develop a method which can simultaneously and iteratively improve data and model performance. A method based on data reconciliation technique has been developed for this purpose. An industrial C3 selective hydrogenation process has been analysed as a case study. The steady state simulator of the process has been developed in Aspen Plus flowsheeting software. Calculated and historical process data were compared using Aspen Simulation Workbook. 74 The proposed technique can be used for monitoring complex processes. Results showed that the proposed technique slightly improves the quality of flow data. Thanks to the accuracy of flow measurements, these process values are directly applicable for the validation of simulator and suitable for further studies, e.g., for the determination of kinetic parameters. Acknowledgements This research of Barbara Farsang and Janos Abonyi was supported by the European Union and the State of Hungary, co-financed by the European Social Fund in the framework of TÁMOP 4.2.4. A/2-11-1-2012-0001 ‘National Excellence Program’. The infrastructure of the research has been supported by the European Union and the Hungarian Republic through the TÁMOP 4.2.2/A- 11/1/KONV-2012-0071 project. SYMBOLS A incidence matrix Ay incidence matrix of balances in terms of measured flows Az incidence matrix of balances in terms of unmeasured flows By Jacobian matrix of measured variables Bz Jacobian matrix of unmeasured variables b source of extensive quantity d y � y̌ error f function vector of equality model constraints f̃ balance error (Aỹ � b vector) f̂ balance error (A˜̂y � b vector) g function vector of inequality model constraints I identity matrix MA methyl-acetylene MAPD methyl-acetylene and propadiene PD propadiene P projection matrix Vd̃ covariance matrix of measurements y̌ vectors of real values of measured variables ỹ vectors of measured variables ŷ vectors of reconciled measured variables ž vectors of real values of unmeasured variables z̃ vectors of unmeasured variables ẑ vectors of estimated unmeasured variables � Lagrange multiplier �Hr reaction heat REFERENCES [1] BAUER, M., CRAIG, I. K. Economic assessment of advanced process control - A survey and framework, Journal of Process Control, 2008, 18(1), 2–18 [2] PACH, F. P. , FEIL, B. , NÉMETH, S., ÁRVA, P. , ABONYI, J. Process data warehousing based operator support system for complex production technologies, IEEE Transactions on Systems, Man and Cybernetics, Part A, Special Issue on ’Advances in Heterogeneous and Complex System Integration’, 2006, 136–153 [3] TIETZE, U., SCHENK, C. Analóg és digitális áramkörök (Műszaki Könyvkiadó) [4] KUHEN, D., DAVIDSON, H. Computer Control II. Mathematics of Control, Chemical Engineering Progress, 1961, 57, 44–47 [5] NARASIMHAN, S. Introduction to data reconciliation and gross error diagnosis (2013), URL: www.che.iitm.ac.in/ naras/ ch544/introDRGED.pdf [6] BELLEC, S., JIANG, T., KERR, B., DIAMOND, M., STUART, P. On-line processing and steady state reconciliation of pulp and paper mill process data, 91st Annual Meeting Preprints-Book B, 2005, 105– 110 [7] THIBAULT, J. Introduction to Data Reconciliation 2013, URL: www.polymtl.ca/namp/docweb/ Modules_Web/M11_Tier1_Chap1-3.pdf [8] NOGITA, S. Statistical Test and Adjustment of Process Data, Ind. Eng. Chem, Process Des. Dev., 1972, 11, 197–200 [9] MAH, R. S. H., STANLEY, G. M., DOWNING, D. M. Reconciliation and Rectification of Process Flow and Inventory Data, Ind. Eng. Chem, Process Des. Dev., 1976, 15, 175–183 [10] MAH, R. S. H. Chemical process structures and information flows, AIChE Journal, 1993, 39(4), 730 [11] CAO, S., RHINEHART, R. An efficient method for on-line identification of steady state, Journal of Process Control, 1995, 5(6), 363–374 [12] TJOA, I. B., BIEGLER, L. T. Simultaneous solution and optimization strategies for parameter estimation of differential � algebraic equation systems, Industrial and Engineering Chemistry Research, 1991, 30, 376–385 [13] BRITT, H. I. , LUECKE, R. H. The estimation of parameters in nonlinear implicit models, Technometrics, 1973, 15, 233–247 [14] MADRON, F., VEVERKA, V. Optimal selection of measuring points in complex plants by linear models, American Institute of Chemical Engineering Journal, 1992, 38(2), 227 [15] CROWE, C. M. Reconciliation of process flowrates by matrix projection. Part II: The nonlinear case, AIChE Journal, 1986, 32, 616–623 75 [16] BAGAJEWICZ, M. J., JIANG, Q. Integral approach to plant linear dynamic reconciliation, AIChE Journal, 1997, 43, 2546–2558 [17] PRATA, D. M., LIMA, E. L., PINTO, J. C. Simultaneous data reconciliation and parameter estimation in bulk polypropylene polymerizations in real time, Macromolecular Symposia, 2006, 243, 91–103 [18] PRATA, D. M., SCHWAAB, M., LIMA, E. L., PINTO, J. C. Nonlinear dynamic data reconciliation and parameter estimation through particle swarm optimization: Application for an industrial polypropylene reactor, Chemical Engineering Science, 2009, 64 18, 3953–3967 [19] MARTÍNEZ-MARADIAGA, M., BRUNO, J. C., CORONAS, A. Steady-state data reconciliation for absorption refrigeration systems, Applied Thermal Engineering, 2013, 51(1-2), 1170–1180 [20] SARABIA, D., DE PRADA, C., GOMEZ, E., GUTIERREZ, G., CRISTEA, S.,SOLA, J., GONZALEZ, R. Data reconciliation and optimal management of hydrogen networks in a petrol refinery, Control Engineering Practice, 2012, 20(4), 343–354 [21] CHEN, J., PENG, Y., MUNOZ, J. C. Correntropy estimator for data reconciliation, Chemical Engineering Science, 2013, 104, 1019–1027 [22] KELLER, T., HOLTBRUEGGE, J., GÓRAK, A. Transesterification of dimethyl carbonate with ethanol in a pilot-scale reactive distillation column, Chemical Engineering Journal, 2012, 180, 309–322 [23] DOMINGUES, L:, PINHEIRO, C. C. , OLIVEIRA, N. M., VILELAS, A., FERNANDES, J., RIBEIRO, F. R. Estimation of catalyst deactivation parameters of ethyl tert-butyl ether (ETBE) reactors based on industrial plant data, in 22nd European Symposium on Computer Aided Process Engineering, Computer Aided Chemical Engineering, vol. 30 (I. D. L. BOGLE, M. FAIRWEATHER, eds.) (Elsevier) 2012, pp. 1002–1006 [24] PRATA, D. M. , SCHWAAB, M., LIMA, E. L., PINTO, J. C. Simultaneous robust data reconciliation and gross error detection through particle swarm optimization for an industrial polypropylene reactor, Chemical Engineering Science, 2010, 65(17), 4943–4954 [25] VASEBI, A., POULIN, E., HODOUIN, D. Dynamic data reconciliation in mineral and metallurgical plants, Annual Reviews in Control, 2012, 36(2), 235–243 [26] VALDETARO, E. D., SCHIRRU, R. Simultaneous Model Selection, Robust Data Reconciliation and Outlier Detection with Swarm Intelligence in a Thermal Reactor Power calculation, Annals of Nuclear Energy, 2011, 38(9), 1820–1832 [27] MARTINS, M. A., AMARO, C. A., SOUZA, L. S., KALID, R. A. , KIPERSTOK, A. New objective function for data reconciliation in water balance from industrial processes, Journal of Cleaner Production, 2010, 18(12), 1184–1189 [28] STANLEY, G., MAH, R. Observability and redundancy in process data estimation, Chemical Engineering Science, 1981, 36(2), 259–272 [29] CROWE, C. M. Observability and redundancy of process data for steady state reconciliation, Chemical Engineering Science, 1989, 44(12), 2909–2917 [30] ALI, Y., NARASIMHAN, S. Redundant sensor network design for linear processes, American Institute of Chemical Engineering Journal, 1995, 41(10), 2237 [31] ALMASY, G.: Mérlegegyenletek - mérési hibák (Egyetemi Kiadó, Veszprém) 1996 [32] NAKAYA, M., LI, X. On-line tracking simulator with a hybrid of physical and Just-In-Time models, Journal of Process Control, 2013, 23(2), 171–178 [33] THIBAULT, J. tate Data reconciliation for Bilinear Systems, 2013, URL: www.polymtl.ca/namp/docweb/ Modules_Web/M11_Tier1_Chap4-5.pdf [34] FLORES, A.: Introduction to Aspen Plus (12. 09. 2001) [35] KANOA, M., NAKAGAWA, Y. Data-based process monitoring, process control and quality improvement: Recent developments and applications in steel industry, Computers and Chemical Engineering, 2008, 32, 12–24 [36] KADLEC, P., GABRYS, B., STRANDT, S. Data- driven Soft Sensors in the process industry, Computers and Chemical Engineering, 2009, 33, 795–814 Advertise upcoming meetings, conferences, and workshops; make public announcements; introduce your research laboratory; a new product or a service in the Hungarian Journal of Industry and Chemistry Please contact us if interested! EDITORIAL OFFICE: UNIVERSITY OF PANNONIA P.O. BOX 158, VESZPRÉM H-8201 (HUNGARY) Tel.: +36 (88) 624-298, E-mail: hjic@almos.uni-pannon.hu; web: hjic.mk.uni-pannon.hu Felelős szerkesztő: Szilágyi Róbert, PhD Kiadja: Pannon Egyetem, 8200 Veszprém, Egyetem u. 10. Levélcím: H-8201 Veszprém, Postafiók 158, Tel.: (88) 624-000 HUNGARIAN JOURNAL OF INDUSTRY AND CHEMISTRY VESZPRÉM Vol. 41(1) pp. 77-82 (2013) STATISTICAL PROCESS CONTROL BASED PERFORMANCE EVALUATION OF ON-LINE ANALYSERS TIBOR KULCSÁR1 AND JÁNOS ABONYIB1 1 Department of Process Engineering, University of Pannonia, Egyetem Street 10., H-8200 Veszprém, HUNGARY BE-mail: janos@abonyilab.com On-line analyzers can provide accurate and timely information for process control and monitoring. Statistical Process Control (SPC) techniques can be effectively utilized to support the development and maintenance of these tools. The D6299- 10 ASTM standard details how on-line analyzers should be validated. The applicability of this standard is demonstrated through the analysis of industrial data collected from an on-line gas chromatograph. The results confirm that automatized SPC can effectively improve the reliability of advanced process control systems. Keywords: on-line analyzer, software sensor, statistical process control, SPC, process monitoring Introduction Process variables characterizing and influencing product quality have a significant role in process control and opti- mization. Off-line laboratory tests mostly take more than two hours. This time delay can cause control problems re- sulting in economic loss. In such situations, an improved on-line monitoring system is required. On-line analyzers eliminate the dependence on laboratory data. Analysers are valuable instruments for real time control because of their fast response time (1-4 minutes) (see Fig. 1) [1]. Quality control techniques can be effectively used to support the development of on-line analyzers [2] and ad- vanced process control systems [3]. The D 6299 ASTM standard (Applying Statistical Quality Assurance Tech- niques to Evaluate Analytical Measurement System Per- formance) provides information for the design and oper- ation statistical quality control (QC) tools to monitor and control of analytical measurement systems using a col- lection of statistical quality control (SQC) tools [4]. The goal of the performance monitoring is the peri- Figure 1: Soft sensors and on-line analyzers enable feedback control odic comparison of the on-line analyzer’s results to the reference value of the same sample measured by labora- tory test methods. Precision and bias (see Fig. 2) are cal- culated to provide information for updating test methods as well as for indicating areas of potential improvements. Control charts and other statistical techniques can be used for performance monitoring. Statistical estimates of the measurement system precision and bias can be calcu- lated on the basis of periodically updated data. Plotting and interpreting these test results can ascertain the in- statistical-control status of the measurement system [5]. On-line Statistical Process Control (SPC) based real-time validation of measurement systems has been already re- Figure 2: Type of faults. The dashed line shows when the fault occurs. �: data free of fault; • corrupted data for the following cases: (a) bias, (b) complete failure, (c) drifting, and (d) precision degradation [4] 78 Figure 3: Illustrative example: MOL aromatic block gas chromatograph IP benzene content (m/m %), Lab, Analyser and APC soft sensor ported in 1997 [6]. The SPC based approach can also be applied for off-line performance evaluation [7]. This ap- proach has a wide range of application areas like steel industry [8], electronic device assembling [9], and build- ings’ energy demand monitoring [10]. The goal of this paper is to present the theoretical background and application details of SPC based perfor- mance evaluation of on-line analyzers. The applicabil- ity of the concepts is illustrated by a case study based on data collected from an on-line gas chromatograph of the MOL Plc (see Fig. 3) . Typical patters that show out of control status of the system are also presented. To check the normality of the residuals, an easily applicable and interpretable tool is proposed. The developed demon- stration tools are available at the website of the authors (www.abonyilab.com). Control Chart based Evaluation of System Performance The studied D6299 practice is devoted to a special test- ing of analyzers [4]. Quality Control (QC) test specimen samples from a specific lot are introduced and tested in the analytical measurement system on a regular basis to establish system performance history in terms of both sta- bility and precision. The control chart is one of the seven basic tools of QC. Control charts - also known as SHEWHART charts - are used to determine if process is in a state of sta- tistical control. The analysis of the control chart indi- cates whether the process is currently under control. This means that these charts are used to check the stability of the production. In stable operation, the variations of the process and quality variables are only random, normally distributed variables. In these cases, no corrections of the control parameters are needed. When the chart indicates that the monitored process is out-of-control, the analysis of the chart can help to determine the sources of the vari- ation. Typically, control charts are used for time-series data, though they can be used for data that have logi- cal comparability [11]. In this section, we present these charts and detail how these should be applied for the per- formance assessment of on-line analyzers. Control Charts Run Chart The run chart is a plot of sample values in chronological order (Fig. 4). The run chart can be used to screen data for unusual patterns such as continuous trending in either direction, unusual clustering, and cycles. The run chart of the data used in our case study is shown in Fig. 3. The plotted time series shows the signal of the on-line ana- lyzer and the related laboratory measurements. I Chart The I (individual) chart is a run chart to which control limits and center line have been added (see top panel of Fig. 5). The center line is based on the mean of the sam- ples, Ī = Pn i=1 Ii n (1) Figure 4: Example for a run chart 79 Figure 5: I and MR charts show out of control samples of the gas chromatograph IP benzene content (m/m) % while the upper and lower control limits are based on the estimated variance (range of the data) MR [12] MR = Pn�1 i=1 |Ii+1 � Ii| n� 1 , (2) UCLI = Ī + 2.66MR, (3) LCLI = Ī � 2.66MR. (4) Individual values that are outside the upper or lower con- trol limits are indications of an unstable system, and ef- forts should be made to determine the cause [5]. Option- ally, any one of the following occurrences should be con- sidered as potential signs of instability: 1. Two out of three consecutive results on the I chart that are more than 1.77MR distant from the center line in the same direction; 2. Five consecutive results on the I chart that are more than 0.89MR distant from the center line in the same direction; 3. Eight or more consecutive points in the I chart that fall on the same side of the center line. MR Chart MR (Moving Range) charts are also used to detect un- usual patterns by plotting the sequential range of two val- ues given by �i = MRi = |Ii � Ii�1| , (5) and connecting each point (see Fig. 6 for an example). There is no lower control limit for an MR chart [5]. The upper control limit for the MR chart is given by UCLMR = 3.27MR. (6) Figure 6: Example for a MR chart Figure 7: Example for an EWMA chart EWMA chart An EWMA (Exponentially Weighted Moving Average) chart is used to enhance the sensitivity in detecting mean shifts that are small relative to the measurement system precision (see Fig. 7 for an example). Each EWMA value is a weighted average of the current result and pre- vious results with the weights decreasing exponentially with the age of the reading [5]: EWMA1 = I1, (7) EWMAi = (1� �)EWMAi�1 + �Ii, (8) where � is the exponential weighting factor. For applica- tion of this practice, a � value of 0.4 is recommended. The control limits for the EWMA chart are calculated using a weight (�) as follows: UCL� = Ī + 2.66MR r � 2� � , (9) LCL� = Ī � 2.66MR r � 2� � . (10) The complexity of multivariate and autocorrelated pro- cesses makes it difficult to use standard control charts. To construct simple and interpretable charts, dimensional reduction could also be used as SIMOGLU and MARTIN have done [13], but in case of autocorrelated data, model- based control charts should be applied as KIM and JITPI- TAKLER have done in their research [14]. Pretreatment of Test Results Assessment, control charting, and evaluation are applied only to appropriately pretreated test results. The purpose 80 of pretreatment is to standardize the control chart scales so as to allow for data from multiple check standards to be compared on the same chart. For QC sample test re- sults, no data pretreatment is typically used since results for different QC samples are generally not plotted on the same chart. In our case, the difference between the measure- ments and their accepted reference values (ARVs) are monitored. ARV serves as an agreed-upon reference for comparison and that is derived based on (1) a theoreti- cal value, based on scientific principles, (2) an assigned value, based on experimental work, or (3) a consensus value, based on collaborative experimental work under the auspices of a scientific or engineering group. I = test result�ARV. (11) Assessment of Initial Results In the initial phase of the application assessment tech- niques are applied to test results collected during the startup phase of or after significant modifications to a measurement system. It is required to perform the fol- lowing assessment after at least 15 pretreated results have become available. The purpose of this assessment is to ensure that these results are suitable for deployment of control charts. Pretreated results should first be visually screened for values that are inconsistent with the remainder of the data set, such as those that could have been caused by tran- scription errors. Those flagged as suspicious should be investigated. Discarding data at this stage must be sup- ported by evidence gathered from the investigation. If af- ter discarding suspicious pretreated results there are less than 15 values remaining, collect additional data and start over. The next step is to examine the pretreated results for non-random patterns such as continuous trending in ei- ther direction, unusual clustering, and cycles. One way to do this is to plot the results on a run chart and examine the plot. If any non-random pattern is detected, investi- gate for and eliminate the root cause(s). Typical Control Chart Patterns In the previous sessions, the main charts used for perfor- mance measurement were presented. In SPC, the Western Electric Rules are the decision rules for detecting “out- of-control” or non-random conditions on control charts. Locations of the observations relative to the control chart control limits (typically at ±3 standard deviations) and centerline indicate whether the process in question should be investigated for assignable causes. The Western Elec- tric Rules were codified by a specially-appointed com- mittee of the manufacturing division of the Western Elec- tric Company and appeared in the first edition of its Sta- tistical Quality Control Handbook in 1956 [15]. Their Figure 8: This process is out of control because a point is either above the UCL or below the LCL. For example, in Fig. 5-B at sequence number 28 there is a unique point below the LCL Figure 9: In this case the system produce 2 out of 3 consecutive points either in or beyond zone A. The process on Fig. 5-B shows this behaviour at sequence number 45 where the values are close to the Upper Control Limit. purpose was to ensure that line workers and engineers in- terpret control charts in a uniform way [5, 12]. The eight standard Western electric rules are: 1. The most recent point plots outside one of the 3- sigma control limits (see Fig. 8). If a point lies outside either of these limits, there is only a 0.3% chance that this was caused by the normal process. 2. Two of the three most recent points plot outside and on the same side as one of the 2-sigma control limits (see Fig. 9). The probability that any point will fall outside the warning limit is only 5%. The chances that two out of three points in a row fall outside the warning limit is only about 1%. 3. Four of the five most recent points plot outside and on the same side as one of the 1-sigma control lim- its. In normal processing, 68% of points fall within one sigma of the mean, and 32% fall outside it. The probability that 4 of 5 points fall outside of one sigma is only about 3%. 4. Eight out of the last eight points plot on the same side of the center line, or target value. Sometimes you see this as 9 out of 9, or 7 out of 7. There is an equal chance that any given point will fall above or below the mean. The chances that a point falls on the same side of the mean as the one before it is one in 81 two. The odds that the next point will also fall on the same side of the mean is one in four. The probability of getting eight points on the same side of the mean is only around 1%. 5. Six points in a row increasing or decreasing. The same logic is used here as for rule 4 above. Some- times this rule is changed to seven points rising or falling. 6. Fifteen points in a row within one sigma. In normal operation, 68% of points will fall within one sigma of the mean. The probability that 15 points in a row will do so, is less than 1%. 7. Fourteen points in a row alternating direction. The chances that the second point is always higher than (or always lower than) the preceding point, for all seven pairs is only about 1%. 8. Eight points in a row outside one sigma. Since 68% of points lie within one sigma of the mean, the prob- ability that eight points in a row fall outside of the one-sigma line is less than 1% (see Fig. 10). Normality Checks Since the control chart and limits prescribed in this prac- tice are based on the assumption that the data behavior is adequately modeled by the normal distribution, it is recommended that a test of this normality assumption be conducted. One way to do this is to use a normal proba- bility plot and the ANDERSON-DARLING Statistic [16]. If the results show obvious deviation from normality, the statistical control charting techniques described are not directly applicable to the measurement system [5]. Quantile-quantile plot (q-q plot) is a graphical tool for comparing two probability distributions by plotting their quantiles against each other. The normality plot of the process can be obtained by comparing the empirical dis- tribution of the data against a standard normal distribution Figure 10: Long runs (8 or more consecutive points) either above or below the centerline. (See Fig. 5-A in range 28-40) Figure 11: Normality check of the studied gas chromatograph. (see Fig. 11). When the results are normally distributed, the plot should be approximately linear. Major deviations from linearity are an indication of non-normal distribu- tions of the differences [17, 16]. Conclusions We showed that statistical process control can be effec- tively used to support the development and maintenance of on-line process analyzers. A case study is presented based on the analysis of data taken from the chemical process industry. The proposed concept has been imple- mented in MATLAB. The results illustrate the applica- bility the developed tools in improving the reliability of advanced process control systems. Acknowledgements This research of Janos Abonyi was supported by the European Union and the State of Hungary, co-financed by the European Social Fund in the framework of TA- MOP 4.2.4. A/2-11-1-2012-0001 National Excellence Program. The infrastructure of the research has been sup- ported by the European Union and the Hungarian Repub- lic through the TAMOP 4.2.2/A- 11/1/KONV-2012-0071. The work of Tibor Kulcsar was supported by the Hungar- ian Scientific Research Fund under Grant K77955. REFERENCES [1] GIDWANI, K. K. , BECKMAN, R. F. Evaluation of refinery control systems, ISA Transactions, 1994, 33(3), 217–225 82 [2] ROBERTO, M. F., DEARING, T. I. , MARTIN, S., MARQUARDT, B. J. Integration of Continu- ous Flow Reactors and Online Raman Spectroscopy for Process Optimization, Journal of Pharmaceuti- cal Innovation, 2012, 7(2), 69–75 [3] SMITH, C. Advanced Process Control: Beyond Sin- gle Loop Control (John Wiley & Sons), 2011 [4] A. INTERNATIONAL: Applying Statistical Qual- ity Assurance Techniques to Evaluate Analytical Measurement System Performance (ASTM Interna- tional), 2010 [5] MONTGOMERY, D. Introduction to statistical qual- ity control (John Wiley & Sons, Inc.), 2009 [6] CONVERSE, J. Real time validation of on-line mea- surements utilizing internal standard introduction, ISA transactions, 1997, 36(2), 151–158 [7] COLLEDANI, M., TOLIO, T. Performance evalua- tion of production systems monitored by statisti- cal process control and off-line inspections, Inter- national Journal of Production Economics, 2009, 120(2), 348–367 [8] KANO, M., NAKAGAWA, Y. Data-based process monitoring, process control, and quality improve- ment: Recent developments and applications in steel industry, Computers & Chemical Engineering, 2008, 32(1-2), 12–24 [9] REIS, M. S., DELGADO, P. A large-scale statistical process control approach for the monitoring of elec- tronic devices assemblage, Computers & Chemical Engineering, 2012, 39, 163–169 [10] BRAGA, L., BRAGA, A. R., BRAGA, C. On the characterization and monitoring of building energy demand using statistical process control methodolo- gies, Energy and Buildings, 2013, 65, 205–219 [11] POOTS, A. J. , WOODCOCK, T. Statistical process control for data without inherent order, BMC Med- ical Informatics & Decision Making, 12, 2012 [12] CHANDRA, M. Statistical quality control (CRC Press LLC), 2001 [13] SIMOGLOU, A., MARTIN, E., MORRIS, A. Mul- tivariate statistical process control of an industrial fluidised-bed reactor, Control Engineering Practice, 2000, 8, 893–909 [14] KIM, S. B. , JITPITAKLERT, W., PARK, S.-K. , HWANG, S.-J. Data mining model-based con- trol charts for multivariate and autocorrelated pro- cesses, Expert Systems with Applications, 2012, 39(2), 2073–2081 [15] THOMAS, D., ET AL. Statistical Quality Control Handbook (Western Electric Co., Inc), 2011 [16] THODE, H. C. Testing For Normality (CRC Press LLC), 2002 [17] HARRIS, T., SEPPALA, C., DESBOROUGH, L. A review of performance monitoring and assessment techniques for univariate and multivariate control systems, Journal of Process Control, 1999, 9(1), 1– 17