A Causal Model of the Sustainable Use of Resources: A Case Study on a Woodworking Process

Abstract

Controlling the life cycle of natural resources, from extraction within the design and the production of products to handling waste, is crucial to green growth and is a part of advancing a resource-efficient, circular economy where everything is fully utilised. One way of using resources more efficiently for a greener economy is to design a production process that takes cost and energy savings into account. From this point of view, the goal of the article is to create a causal description of sustainable woodworking—especially using renewable and non-renewable resources—in relation to changes in the concentration levels of CO₂ in the atmosphere. After estimating the partial parameters, this model can be used to predict or simulate different CO₂ concentration levels in the atmosphere—for example, based on the ratio of renewable to non-renewable sources. After a theoretical description, the subsequent practical goal is to identify the optimal settings of wood-milling process parameters for either minimising energy consumption per workpiece and unit variable costs or for maximising the overall customer benefit. For this purpose, a complete factorial design was used, and based on this, the consumption energy (direct cost) optimisation of the production process was supplemented by a profitable production calculation. The effect of reducing variability was verified using a statistical F-test. The impact of minimising energy consumption (economically expressed as the mean profit) was then validated using a Student’s t-test.

1. Introduction

The way society currently produces energy is not sustainable [1]. With the increase in greenhouse gas emissions, the world is witnessing the continuous melting of polar ice caps and frequent extreme weather events due to global warming [2]. These effects of global warming have seriously threatened the lives and productivity of humans; thus, the control and reduction of greenhouse gas emissions has become an urgent problem [2,3]. One of the objectives of the Kyoto Protocol, which has been followed by the European Council (European Green Deal, 2020), is to decarbonise production and reduce energy consumption by moving to renewable energy and improving energy efficiency. According to [4,5], a certain degree of global warming is already inevitable. Thus, warming will continue for decades to come, even if we stabilise the carbon dioxide concentration level in the atmosphere at the current level. However, these concentration levels are now significantly higher than they have been at any point in recent history. Around the beginning of the Industrial Revolution, carbon dioxide (CO₂) concentrations fluctuated at 280 ppm (number of CO₂ molecules per million air molecules) [6,7]. We must avoid increasing concentrations above 450 ppm, and according to the Paris Climate Agreement (2015), it is necessary to limit the temperature increase to 1.5 °C to avoid drastic social changes.

Many scientists [8,9,10,11,12] have noted that a suitable long-term target is to decrease the CO₂ concentration levels to 350 ppm. However, this would require the removal of CO₂ from the atmosphere, whereas the world is currently adding CO₂. The sustainability of the green circular economy will only emerge when society significantly reduces its energy consumption. Efficiency policies are available, but they still do not achieve what is theoretically possible and within reason. Thus, an essential area for improvement is energy efficiency [13]. Studies show that it is possible to reduce the amount of energy consumed by industrial electric sources by 120 TWh per year within the European Union by 2025, which would be enough of a reduction that several central power stations would no longer be needed.

Unfortunately, not enough has been done to boost energy efficiency. Although the European Union has obligatory carbon emission objectives, the energy efficiency target is only voluntary. The only obligatory part is a reduction of greenhouse gas emissions by 40% until 2030. The aim of increasing the share of renewable energy generation for the same period is 27%, but it is only binding for the European Union. There are no specific targets for EU member states. Here, the efficiency target is 27%, but this is also optional. According to [14,15], a certain degree of global warming is already inevitable.

Research Questions

A factorial design can be used to reduce the variability of woodworking, as well as to significantly reduce energy consumption (or increase the profit margin).

To verify the first part of the research question, we create a zero hypothesis:

Hypothesis 1 (H1). 

The full-factorial design of technological factors has no statistically significant effect on production variability.

Hypothesis 1a (H1a). 

The full-factorial design of technological factors has a statistically significant influence on production variability.

To verify the second part of the question, we create a zero hypothesis:

Hypothesis 2 (H2). 

The full-factorial design of technological factors has no statistically significant effect on the mean value of energy consumption.

Hypothesis 2a (H2a). 

The full-factorial design of technological factors has a statistically significant effect on the mean value of energy consumption.

2. Literature Review

Because manufacturing industries are facing energy sustainability challenges because of increasing global competition, they need to continually increase productivity while reducing manufacturing costs. The wood processing industry is a significant consumer of energy and other resources, causing a severe environmental impact [16,17,18]. Hence, decreasing the energy requirement of manufactured products can be a suitable target for enhancing both economic competitiveness and ecological sustainability. Various models have tackled the need to decrease manufacturing energy consumption. ISO 50001, for example, implements a systematic approach to continuously improving energy performance and defining the specifications for process design and documentation [19]. Furthermore, EN 16231 norms recommend a methodology for energy data evaluation to discover the energy effectiveness of particular parts (such as wood-milling devices), enabling electric energy consumption monitoring and correlating this monitoring with other units [20]. Here, energy capability benchmarking is a system that shows energy consumption.

A negative side effect of a higher range of production is that machine tools can be substantial industrial energy consumers [21,22]. Therefore, a reduction in the demand for energy from machine tools can significantly improve the production processes’ environmental impact and the carbon footprint of consumer products. When looking at the carbon footprint, a typical product life cycle consists of three stages: production, application and end of life. The use phase is the most energy-intensive phase for a machine tool itself, causing 55% to 85% of CO₂ emissions during its life cycle. Recently, a design for the environmental evaluation of machine tools during their use phase was introduced. The methodology represents a reproducible quantification of energy consumption in different process settings. Gontarz et al. presented a modular arrangement procedure for machine tools based on multichannel analyses to improve energy efficiency and enable total cost of ownership (TCO) calculations [21].

Several studies have also been carried out to model the energy consumption of machine tools and, thus, determine the environmental impact of the goods produced [23,24,25]. The machining time is a crucial aspect of the energy consumption of machine tools, especially those with a high baseload (i.e., machines with many supporting parts such as hydraulic, exhaust and cooling lubricant systems, etc.). Neugebauer et al. focused on system-level events for energy-efficient increase of machine tools and creating a production system with a direct efficiency increase on the component level using an optimised interface of the components on the higher system level [5]. Mori et al. emphasised that an energy consumption reduction can be obtained by modifying cutting conditions for regular drilling, face/end milling and various machining operations [25]. The overall influence of lightweight design methods on the energy efficiency of machine tools and limitations on the maximum mass reduction for structural components and lightweight machine tools was studied by Kroll et al. [26].

Various researchers have confirmed that high material removal rates generally reduce machine tools; energy consumption because of the reducing machining time [17,25,27]. A typical energy performance indicator (EnPI) used for benchmarking within or between units is the specific energy consumption (such as the energy per unit produced). Emerging trends in these indicators help us to validate changes in energy efficiency, but also act as evidence for issues such as process plan deviations and changes in process stability or quality.

Further optimisation potential arises from a proper choice of the tool path strategy during machining [8,24,27]. As a conclusion, it is of high importance to use optimal machining procedures and parameters in combination with performant tooling systems in order to minimise the cycle times and thus the energy consumption. However, over the course of machining, the process performance might change due to tool wear, suboptimal machine settings or operating errors. Hence, it is expedient to monitor certain performance indicators over time in order to assess and compare different processes [1]. It is important to use optimal machining procedures and parameters in combination with performance tooling systems to minimise cycle times and, thus, energy consumption. However, over the course of machining, the process performance might change because of tool wear, suboptimal machine settings or operating errors. Hence, it is crucial to be able to monitor certain performance indicators over time to assess and compare different processes.

Principle of Milling

Milling is a machining operation in which a layer of material, in the form of small, individual chips, is taken from a workpiece by a rotary multitooth tool—a milling cutter. The milling cutter rotates around its axis while working, and gradually cuts into a workpiece, which is simultaneously moved against the tool. Each cutter tooth gradually cuts short chips of uneven thickness from the material to be machined so that the cutting process is intermittent. Using this method and various types of milling tools, it is possible to machine on workpieces, mainly planar surfaces, but also shaped, oblique, irregular and rotary surfaces, dividing the material into different lengths. This wide application and the possibility of precise production have given milling an important role in engineering. In most cases, high-speed milling allows more productive and cost-effective material removal than single-edge machining such as turning or planning. In some complex cases, milling is the only machining that can be used [28].

The electrical energy efficiency of the deformation work during the main cutting movement is 45% (±20%) [29], which is higher than the mean value of 32% (±15%) for single cutting. An important advantage of milling is the higher shape matching over single-point machining—for example, turning [30,31]. Several methods can be used to measure the energy machining efficiency (EME). The specific cutting energy (SCE) is described in [32] as the energy spent per unit of material removal volume (J/cm³). The amount of specific energy is determined by the material’s ability to form chips, blade geometry and material hardness. The EME has previously been investigated, for example, using experimental methods of analysis [33,34], and it has been studied [35] for the selection of coolant using a design of experiment–response surface methods. A more specific experiment applied to steels of different alloy and carbon contents was performed in [15]. Lee et al. [24] presented another research direction for the study of milling factors that affect the residual stress and surface roughness. Some later investigations focussed on the effect of EME and SCE on undeformed chip thicknesses (i.e., machining without residual stress strain). Altintas [36] designed an EME/SCE regression model for milling using a Taguchi experiment. Methodically, other EME/SCE research has been presented [37,38] based on a fuzzy logic approach developed to optimise the cutting parameters in micromilling.

The abovementioned studies are sensitive to the proper design of the experiment—that is, using certain principles to eliminate random and systematic error (randomisation, blocking, replication, design balance, etc.). An occasional inability to fulfil these principles slightly limits the application of the proposed methodologies. In this respect, EME/SCE analytical models are coming to the forefront. The development of analytical methods has shifted from the tool life phase to the current phase of factor optimisation. In addition to empirical models, analytical methods can also be used for energy-saving milling. For example, Pawade et al. [33] and Wang [31] designed models to specify the effect of cutting parameters on EME/SCE. Han [13] created a logical milling model to optimise surface roughness. In recent years, it has proven to be energy-efficient to use vortex milling [34], which has approximately nine times the cutting rate of conventional milling, as well as a 25% lower cost and comparable surface quality [13,39,40].

3. Materials and Methods

The planning and operation of the energy system requires a sustainability assessment of the system, in which the load model adopted is the most important factor [34]. In the woodworking industry, it can be useful to consider sustainable development as a system and to investigate its related behaviour by looking at the solution to the system task.

3.1. System Concept of the Functioning of Sustainable Development in Terms of CO₂ Balance in the Atmosphere

It is evident from Figure 1 that the woodworking industry fulfils its goals through the wooden products offered. This is the basis of the general transformation process, which is implemented with the participation of all available renewable resources: energy (e_r) material (h_r), transferable knowledge (k_r) and personnel (p_r). Renewable energy resources are defined as ‘energy generated from natural resources that can be naturally replenished in the environment’ through sustainable energy resources. These resources include hydroelectric, wind, biomass, geothermal and solar energy.

The system output is the response to the stimulus of the input impulses from the non-renewable environment (non-renewable resources) and is composed of matter (h_n), fuel energies (e_n) and tacit (non-transferable) knowledge (k_n). The substances represent the material and raw materials transformed into the form of the requested goods; energy is used to implement the transformation process, and knowledge shapes the transformation process so that the outputs are competitive. The open character of the woodworking industry system is achieved by its interaction with the environment, thus, the input/output behaviour of the system is causally unstable.

Resource consumption management is implemented by an imaginary resource pump that uses the distributor to push the renewable and non-renewable resources into the transformation process, which implements the conversion of inputs into an output—the product.

3.2. Input–Output Model in the Case of a Stable Condition

In a stable condition, the amount of emitted CO₂ (and possibly other greenhouse gases, including CH₄ (methane), N₂O (nitrous oxide) and CFC (freons)) is fully absorbed by the oceans and absorbed, for example, by liquid solvents. Thus, in the steady state, there is no increase in CO₂ in the atmosphere (ΔCO₂ = 0). In this case, the change in the value of CO₂ can be expressed by the total differential, consisting of partial differential equations (see Equation (5)).

In accordance with Figure 1, the total change in the CO₂ emitted at the unit consumption of a given renewable or non-renewable source is the sum of their partial changes. For a relatively short time when there is an increase in the i-type resources consumed and Δqi only has low change values, it is possible to linearise the dependent variable function (increase in CO₂ emissions) based on a relatively small time change. This will cause only a small deviation from the real functional value.

Now, it is possible to express the change in the value of CO₂ emissions for the time interval ∆t as the direction of the function CO₂ = f(t) at point CO₂ (derivation of the function in the point), multiplied by the time interval Δt, that is:

$${\mathsf{\Delta }\mathrm{CO}}_2=\frac{\mathrm{d}\left({\mathrm{CO}}_2\right)}{\mathrm{dt}}\times \mathsf{\Delta }\mathrm{t}.$$

(1)

If expressing this change in the value of CO₂ emissions depends on the amount of resources consumed in a certain time, the linearisation of the incremental value of ΔCO₂ will be determined by the following equation:

$${\mathsf{\Delta }\mathrm{CO}}_2=\frac{\mathrm{d}\left({\mathrm{CO}}_2\right)}{{\mathrm{dq}}_{\mathrm{e}}}\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{e}}.$$

(2)

In the same way, further increases in CO₂ emission values could be linearised depending on the volume of the supplied resource during a certain period. Then, equations can be created whereby it is possible to describe the behaviour of the wood processing industry. The greatest increases in CO₂ are represented by non-renewable material and energy sources. However, even renewable sources represent CO₂ emissions (e.g., in the production of equipment—a power plant powered by a renewable source). Even the creation of tacit knowledge is a marginal source of CO₂ growth [40]. For example, to create knowledge about the technological process of wood processing, a set of experiments must be performed in which matter is transformed using energy.

Expressing the change in ∆CO_2n emissions caused by the consumption of non-renewable resources is equal to the sum of partial changes for individual types of non-renewable resources that are used. This concerns changes in matter (hn), fuel energies (en) and tacit (non-transferable) knowledge (kn) resources. When linearising individual non-renewable resources, it is possible to express the change in ∆CO_2n emissions caused by the total differential:

$${\mathsf{\Delta }\mathrm{CO}}_{2\mathrm{n}}=\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{en}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{en}}+\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{hn}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{hn}}+\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{kn}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{kn}}.$$

(3)

When linearising individuals by the consumption of renewable resources, it is possible to characterise the change in the value ∆CO_2r emissions caused by the consumption of renewable resources by Equation (4):

$${\mathsf{\Delta }\mathrm{CO}}_{2\mathrm{r}}=\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{er}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{er}}+\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{hr}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{hr}}+\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{kr}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{kr} }+ \left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{pr}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{pr}}.$$

(4)

During stable conditions, a change in the ∆CO₂(ABT) absorption is balanced by the increase in CO₂ from the use of renewable (∆CO_2r) and non-renewable (∆CO_2n) sources:

$$\Delta \mathrm{CO}2\left(\mathrm{ABT}\right) =\Delta \mathrm{CO}2\mathrm{r} +\Delta \mathrm{CO}2\mathrm{n},$$

(5)

The stable condition means a zero increase in CO₂ emissions into the atmosphere (change in ΔCO₂ = 0, so that production in the woodworking industry is carbon-neutral). The individual variables ∆CO_2r and ∆CO_2n are replaced by their linearised components, stated by Equations (3)–(5), resulting in:

$$\begin{gathered} {\left(\frac{\partial \mathrm{CO}2{\left(\mathrm{ABT}\right)}_{}}{\partial {\mathrm{q}}_{\left(\mathrm{ABT}\right)}}\right)}_0\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{ABT}}=\left[\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{en}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{en}}+\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{hn}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{hn}}+\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{kn}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{kn}}\right]+ \\ +\left[\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{er}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{er}}+\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{hr}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{hr}}+\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{kr}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{kr} }+ \left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{pr}}}\right)\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{pr}}\right]. \end{gathered}$$

(6)

When modifying Equation (6) so that the individual addend from both square brackets is merged on the right side of the equation, it results in Equation (7), in which there is a generic grouping of individual sources:

$$\begin{gathered} {\left(\frac{\partial \mathrm{CO}2{\left(\mathrm{ABT}\right)}_{}}{\partial {\mathrm{q}}_{\left(\mathrm{ABT}\right)}}\right)}_0\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{ABT}}={\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{en}}}+\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{er}}}\right)}_0\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{e}}+{\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{hn}}}+\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{hr}}}\right)}_0\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{h}}+ \\ +{\left(\frac{\partial {\mathrm{CO}}_2}{\partial \mathrm{kn}}+\frac{\partial {\mathrm{CO}}_2}{\partial \mathrm{kr}}\right)}_0\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{k}}+{\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{pr}}}\right)}_0\times {\mathsf{\Delta }\mathrm{q}}_{\mathrm{p}}. \end{gathered}$$

(7)

From Equation (7), the importance of the efficient use of energy sources in the wood processing industry is shown via the important and replaceable role of renewable resources within the total source relations of the industry. The optimising (minimising) of the partial derivation of the value of the non-renewable consumption according to the volume of this source $\left(\frac{\partial {\mathrm{CO}}_2}{\partial {\mathrm{q}}_{\mathrm{n}}}\right)=\mathrm{Minimum}$ is a key issue in cases where we cannot totally replace non-renewable resources.

Leaving behind the idealised situation of the CO₂ levels in the atmosphere, where the behaviour of the wood processing system in a stable condition was derived, in a real situation, there is an increase in CO₂ levels in the atmosphere because not all of the emitted amounts are absorbed. ∆CO₂(ABT) absorption is balanced with the increase in CO₂ from the use of renewable (∆CO_2r) and non-renewable (∆CO_2n) sources. Therefore, the value of ∆CO₂(ABT) absorption is not equal to the increase in CO₂ levels from the use of renewable (∆CO_2r) and non-renewable (∆CO_2n) sources. In this case, it is possible to express this value difference as follows:

$$\Delta \mathrm{CO}2\mathrm{r} +\Delta \mathrm{CO}2\mathrm{n}- \Delta \mathrm{CO}2\left(\mathrm{ABT}\right)=\frac{\partial \left(\mathrm{CO}2\mathsf{\Delta }\mathrm{hz}\right)}{\partial \mathrm{t}}\times \left({\mathsf{\Delta }\mathrm{q}}_{\mathrm{r}}+{\mathsf{\Delta }\mathrm{q}}_{\mathrm{n}}\right),$$

(8)

where the expression $\frac{\partial \left(\mathrm{CO}2\mathsf{\Delta }\mathrm{hz}\right)}{\partial \mathrm{t}}$ is the immediate change in the accumulated CO₂ emissions at a certain point in time (the year 2019 already reached values higher than 400 ppm).

Equation (8) is a first-order differential equation. From the results, we get the predictive equation of CO₂ concentrations in the atmosphere, which is based on a simple causal description. Unlike the extrapolation forecast, it does not need to know retrospective data to predict future developments.

3.3. Energy Optimisation of Milling (with the Economic Limitation of Profitability of Production)

From the start of civilisation, human beings have been involved in the process of transforming natural resources into useful products. There are many machining operations that transform raw materials into useful products. Milling is efficient as well as flexible, and thus, has been widely adopted by most manufacturing industries. These manufacturing industries are responsible for the most conversion of natural resources. Hence, a more sophisticated way of optimising the manufacturing process is required and the appropriate selection of various operating parameters is necessary [2]. Minimum energy consumption, desirable surface quality and maximum material removal are some of the common challenges faced by all types of industries and influence the cost of production as well as the quality of products [5,35,40]. Hence, a more methodical approach is required that uses experimental methods, mathematical models and statistical analysis to bring down the cost of production by optimising these operating parameters [38,39].

The standard parameter of the machine is its maximum power input, which can be expressed in kVA (apparent power), according to which circuit breakers affect the electricity fees. In practice, machine components never run at maximum power simultaneously. Thus, if we sum up their performance, we will obtain an unrealistically high number. This can be illustrated by servo drives interpolating axes that reach peak power only for a short time and at different moments than the spindle does [33]. Therefore, the value of the maximum power input can be empirically reduced to represent a limit that will not be exceeded in practice [13,41]. The installed power is still far from the mean power input—that is, the arithmetic mean over a longer working time. This estimation uses either the nominal working power input of the machine components or its distribution over time. The medium outputs determine the importance of individual appliances and their priority in energy optimisation. Sometimes, for example, the extraction of fog from the working area becomes more important than the drivers of all the linear axes (e.g., in a small milling centre), even though it has a fractional installed capacity. For a correct determination of the mean wattages of machine components, it is best to perform measurements with a multichannel wattmeter during typical operating modes.

It is possible to optimise the energy and cost production process in terms of setting the production parameters [42]. This optimisation is based on contradictory factors: a more energy-efficient mode of production is often accompanied by a lower efficiency of time capacity and, consequently, the lower productivity of the production equipment. This case leads to an increase in the average fixed costs. Therefore, in practice, it is a question of finding an optimal ratio between the average variable costs and average fixed costs. At a constant work price of the workpiece (i.e., average yield), this finding leads to maximising the profit margin.

The progress includes the development of a full-factorial model that covers the relationship into the empirical mean profit (per workpiece) and yield from the experimental trial, which is determined by significant woodworking process parameters. This factorial model is based on the following statistical relevant factors: workpiece speed (nw), cutting speed (nt) and axial feed speed (vf). The estimated ratio of the mean profit of identical products $\widehat{y}$ is given by the sequential equation of a regression model for three factors [43]:

$$\hat{\mathrm{y}}={ \mathsf{\beta }}_0+{\mathsf{\beta }}_1{\mathrm{x}}_1+{\mathsf{\beta }}_2{\mathrm{x}}_2+\cdots +{ \mathsf{\beta }}_{12}{\mathrm{x}}_1{\mathrm{x}}_2+{\mathsf{\beta }}_{13}{\mathrm{x}}_1{\mathrm{x}}_3+\cdots +\epsilon ,$$

(9)

where ${\mathsf{\beta }}_1$ and ${\mathsf{\beta }}_2$ are the regression coefficients (calculated as half the effect of a given factor) and ${\mathsf{\beta }}_0$ is the distance of the origin of the response surface at the beginning of the Cartesian system (calculated as the average response in the factorial experiment). The term $\mathsf{\epsilon }$ is the normally distributed random error. The regression coefficients ${\mathsf{\beta }}_{12}$ and ${\mathsf{\beta }}_{13}$ correspond to the interaction between the process parameters x₁ and x₂ and x₁ and x₃, respectively.

If we label the three factors of the two-level, full-factorial design as A, B and C, then the effect of each factor and interactions between the factors can be formally calculated according to the following equations [42]:

Estimation of the effect of factor A:

$$\mathrm{A}={\overline{\mathrm{y}}}_{{\mathrm{A}}^{+}}-{\overline{\mathrm{y}}}_{{\mathrm{A}}^{-}}=\frac{1}{4\mathrm{n}}\left[\mathrm{a}+\mathrm{ab}+\mathrm{ac}+\mathrm{abc}-\mathrm{b}-\mathrm{c}-\mathrm{bc}-\left(1\right)\right].$$

(10)

Estimation of the effect of factor B:

$$\mathrm{B}={\overline{\mathrm{y}}}_{{\mathrm{B}}^{+}}-{\overline{\mathrm{y}}}_{{\mathrm{B}}^{-}}=\frac{1}{4\mathrm{n}}\left[\mathrm{b}+\mathrm{ab}+\mathrm{bc}+\mathrm{abc}-\mathrm{a}-\mathrm{c}-\mathrm{ac}-\left(1\right)\right].$$

(11)

Estimation of the effect of factor C:

$$\mathrm{C}={\overline{\mathrm{y}}}_{{\mathrm{C}}^{+}}-{\overline{\mathrm{y}}}_{{\mathrm{C}}^{-}}=\frac{1}{4\mathrm{n}}\left[\mathrm{c}+\mathrm{ac}+\mathrm{bc}+\mathrm{abc}-\mathrm{a}-\mathrm{b}-\mathrm{ab}-\left(1\right)\right].$$

(12)

Estimation of the effect of interaction between factors A and B:

$$\mathrm{AB}={\overline{\mathrm{y}}}_{{\mathrm{AB}}^{+}}-{\overline{\mathrm{y}}}_{{\mathrm{AB}}^{-}}=\frac{1}{4\mathrm{n}}\left[\mathrm{ab}+\left(1\right)+\mathrm{abc}+\mathrm{c}-\mathrm{b}-\mathrm{a}-\mathrm{bc}-\mathrm{ac}\right].$$

(13)

Estimation of the effect of the interaction between factors A and C:

$$\mathrm{AC}={\overline{\mathrm{y}}}_{{\mathrm{AC}}^{+}}-{\overline{\mathrm{y}}}_{{\mathrm{AC}}^{-}}=\frac{1}{4\mathrm{n}}\left[\mathrm{ac}+\left(1\right)+\mathrm{abc}+\mathrm{b}-\mathrm{a}-\mathrm{c}-\mathrm{ab}-\mathrm{bc}\right].$$

(14)

Estimation of the effect of the interaction between factors B and C:

$$\mathrm{BC}={\overline{\mathrm{y}}}_{{\mathrm{BC}}^{+}}-{\overline{\mathrm{y}}}_{{\mathrm{BC}}^{-}}=\frac{1}{4\mathrm{n}}\left[\mathrm{bc}+\left(1\right)+\mathrm{abc}+\mathrm{a}-\mathrm{b}-\mathrm{c}-\mathrm{ab}-\mathrm{ac}\right].$$

(15)

Estimation of the effect of the interaction between factors A, B and C:

$$\mathrm{ABC}={\overline{\mathrm{y}}}_{{\mathrm{ABC}}^{+}}-{\overline{\mathrm{y}}}_{{\mathrm{ABC}}^{-}}=\frac{1}{4\mathrm{n}}\left[\mathrm{abc}-\mathrm{bc}-\mathrm{ac}+\mathrm{c}-\mathrm{ab}+\mathrm{b}+\mathrm{a}-\left(1\right)\right],$$

(16)

where [42]: ${\overline{y}}_{A^{+}}$ = the mean response factor for the upper-level A, the average response for the lower level of factor A; ${\overline{y}}_{A^{-}}$ = the mean response factor for the lower-level A, the average response for the lower level of factor A; (1), a, b and c = all eight combinations of the responses for two setting levels of the three factors, and n = the number of replications of this design.

In the machining facility where the experiments were performed, a high-speed CNC milling machine SKS-7 GS 2215 HS, Prague [CZ], was used (see Figure 2).

The active electric power P of the alternating electric current (in

Table 1. Measurement of the value of electric energy consumption for the two milling tools.

Electrical Variable	Unit	Average Response Tool: Glocken Messer ADE-06	Average Response Tool: Milling Cutter MKS-V25
Phase current (IL)	kWsec	2.59	2.42
Mechanical power of the engine (PM)	W	1646	1767
Electrical input power (PE)	W	2352	2291
Energy consumed per four products produced (Wc)	kWsec	2893	1710
Time	s	1150	750
Spindle torque (M)	Nm	65.8	62.4
Revolution (f)	min−1	1500	1700
Energy efficiency of power input (η)	%	69.9	77.1

, this is replaced by mechanical power (PE)) is determined by Equation (17):

$$\mathrm{P}= \frac{1}{\mathrm{T}}{{\displaystyle \int }}_0^{\mathrm{T}}\mathrm{p}\times \mathrm{dt}=\frac{1}{\mathrm{T}}{{\displaystyle \int }}_0^{\mathrm{T}}\mathrm{u}\times \mathrm{i} \mathrm{dt},$$

(17)

where p is the instantaneous power (for instance, the instantaneous electric voltage) and the instantaneous electric current and T is the measurement period. The electric active power characterises the irreversible conversion of energy into useful energy (mechanical work of milling).

An electric motor of the milling machine is connected as a load (L) to perform the main milling movement. The Gossen A2000 multifunction wattmeter, Prague [CZ], (W) was used to measure and record the instantaneous power. The wattmeter measures the instantaneous voltage and current values of three phases (R, S and T). In addition, a VoltCraft Vc-610Bt voltmeter (A), Prague [CZ], was used in the role of an ammeter, and a GW Instek GDM-8341 voltmeter (V), Prague [CZ], was used to correct the voltage drop because of wattage consumption. On the given production equipment (Copy Milling Machine COSMA-SKS 7), we performed an electricity consumption per production batch experiment (one production cycle of four workpieces). We also determined the mechanical power of the engine for determining the efficiency of the production machine (η). We performed the reaction torque measurement with a DIS CV-505 dynamometer (DM), Prague, [CZ]. For this characteristic, we measured the reaction moment. When a force is applied at the end of the arm, the beam is deformed and detected by the strain gauge (its electrical resistance varies depending on the magnitude of the deformation). The changed resistance is then converted to voltage, and the torque value is shown directly on the display and stored in memory. The mean torque value and the corresponding spindle speed then determine the motor’s mechanical power (P_M) according to Equation (18), which was further used to assess the efficiency of the milling process (see Figure 3 and Table 1):

$${\mathrm{P}}_{\mathrm{M}}=\frac{\mathrm{M}\times \mathrm{f}}{60},$$

(18)

where P_M is the mechanical power of the engine (W Watt), M is the spindle torque (Nm Newton meter) and f is the steady revolution per minute.

The measured values of the electrical quantities are represented in Table 1. The production cycle’s energy consumption was determined for two types of tools (Glocken Messer ADE-06, Prague [CZ], and the other unique tool developed by the company where the experiment was performed—MKS-V25 milling cutter, Prague [CZ]). Different settings of the production parameters were made for each tool. After comparing the values, the production cycle was significantly more energy-efficient with the MKS-V25 milling cutter tool. Thus, innovations in the form of the development of a specialised milling tool have proven to be meaningful.

For the MKS-V25, we subsequently performed an optimised design of the process parameters in terms of minimising machine time in connection with reliability and productivity—that is, according to the criterion of economic profitability. Specifically, a factor optimisation of three response characteristics were measured using the following scale: productivity: q (has a positive effect on machinery investment); reliability: r (measured as the proportion of compliant production to total production at a given set of process parameters); mean variable costs: vc (measured as the sum of energy costs and material costs per unit of product, see

Table 2. Setting of the variables in absolute and coded values.

Process Parameter	Unit	Low Setting	High Setting	Low Setting (Coded Units)	High Setting (Coded Units)
Axial speed (vf)	mmpsec	0.5	2.0	−1	+1
Workpiece speed (nw)	rpm	2	8	−1	+1
Cutting speed (nt)	rpm	150	200	−1	+1

).

The first purpose of the factorial design was to find the factors and interactions influencing the mean variability of the ratio of the identical products. The result of these trials is shown in Table 2. For the significance test, it was decided to use the level of importance a = 5% (0.05). Then, if the p-value was less than the level of importance (0.05), the factor or interaction would be statistically significant. Productivity (q) has the main effects if the axial speed (vf), workpiece speed (nw), cutting speed (nt) and interactions between these process parameters are statistically significant. This finding is further supported by a normal plot (see Figure 2). Reliability (r) has the main effects if the axial speed (vf) and workpiece speed (nw) are statistically significant (see Figure 3). The mean variable cost (vc) has the main effects if the axial speed (vf), workpiece speed (nw) and cutting speed (nt) process parameters are statistically significant (see Figure 4). In addition, the mean profit (mean TP) shows no statistical dependence on any process parameter or interaction (see Figure 5), and the process parameters are statistically significant. This is probably because of the mutual compensation of production reliability and productivity. Therefore, we determined the mean profit (mean TP) for each process factor setting trial according to the following equation:

$$\begin{gathered} \mathrm{mean} \left(\mathrm{TP}\right)={\mathrm{r}}_{\mathrm{i}}\times \mathrm{p}-\left({\mathrm{vc}}_{\mathrm{i}}+\raisebox{1ex}{$\mathrm{FC}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{q}}_{\mathrm{i}}$}\right.\right), \\ \mathrm{i}\in \mathrm{k}\times \left\{1,2,3,4,5,6,7,8\right\} \end{gathered}$$

(19)

where r_1,2 is the reliability in achieving a suitable workpiece with a given set of process parameters, vc_i is the mean variable costs with a given set of process parameters, p is the mean workpiece price; FC is the fixed costs (expressed as investment in production equipment), q_i = the productivity with a given set of process parameters and k is the number of replications of each combination of process parameter settings (in our case, k = 2).

4. Results

We performed 16 simulation trials based on the two-level, full-factorial design. The normal plots from the design are shown in Figure 4, Figure 5, Figure 6 and Figure 7. The normal probability plot of the effects shows the standardised effects relative to a distribution fit line for the case when all the effects are zero [31]. The standardised effects are t-statistics that test the null hypothesis that the effect is 0. Positive effects increase the response when the settings change from the low value of the factor to the high value. Negative effects decrease the response when they settings change from the low value of the factor to the high value of the factor. Effects further from 0 on the x-axis have greater magnitude. Effects further from 0 are more statistically significant.

The distance that points must be from the reference line to be statistically significant depends on the significance level (denoted by α = 0.05). Unless you use a stepwise selection method, the significance level is 1 minus the confidence level for the analysis.

We use a normal probability plot of the effects to determine the magnitude, direction and the importance of the effects (see Figure 4, Figure 5, Figure 6 and Figure 7). On the normal probability plot of the effects, effects that are further from zero are statistically significant. The colour and shape of the points differ between statistically significant and statistically insignificant effects. For example, in Figure 4, the main effects for factors A (working speed), B (cutting speed) and C (axial speed) and interactions AB, AC and BC are statistically significant at the 0.05 level (the ABC interaction is not statistically significant). These points have a different colour (compared to the ABC interaction) and shape from the points for the insignificant effects.

In addition, the plot indicates the direction of the effect. Processes A (working speed), B (cutting speed) and C (axial speed), and interactions AC, AB and BC, have a positive standardised effect. When these processes change from the low level to the high level of the factor, the response increases. The interaction ABC (between working speed, cutting speed and axial speed) (B) has negative standardised effects (but is not statistically significant, thus, it can be excluded from the optimised process). When the ABC interaction increases, the q (productivity) response decreases.

In other worlds, the terms with the highest positive effect are on the right. According to Figure 4 (productivity response—q), the effects came from the two-way and three-way interaction terms between the ABC factors (axial speed‒workpiece speed‒cutting speed), AB factors (axial speed‒workpiece speed) and BC factors (working speed‒cutting speed). According to Figure 4, the separate factors are also crucial for productivity q. The two-way and three-way interaction terms are not important for further responses (reliability r, see Figure 5; variable cost vc, see Figure 6; mean profit, see Figure 7). Therefore, it is appropriate to use multiple linear models to estimate the mean profit and reliability responses. The interaction causes a ‘curvature’ of the response space, and this means it would be appropriate to use a higher degree model (e.g., polynomial regression).

The calculated effect factors in the coded values (response factor to change from −1 to +1) are listed in the second column in

Table 3. Estimated effects and coefficients for mean profit (coded units).

Term	Effect	Coef.	SE Coef.	t-Value	p-Value	VIF
Constant		3.064	0.118	26.02	0.024
nw (rpm)	−0.155	−0.078	0.118	−0.66	0.629	1.00
nt (rpm)	−0.038	−0.019	0.118	-0.16	0.897	1.00
vf (mmprev)	−0.094	−0.047	0.118	−0.40	0.758	1.00
nw (rpm) × nt (rpm)	0.154	0.077	0.118	0.65	0.631	1.00
nw (rpm) × vf (mmprev)	0.195	0.098	0.118	0.83	0.559	1.00
nt (rpm) × vf (mmprev)	0.175	0.088	0.118	0.74	0.593	1.00
Model Summary
S				R-sq
0.33304				79.59%
Regression Equation in Uncoded Units
Mean profit = 3038 − 0.1287 nw + 0.007000 nt − 0.09825 vf + 0.1283 nw × nt + 0.1715 nw × vf + 0.1387 nt × vf − 0.1915 nw × nt × vf

,

Table 4. Estimated effects and coefficients for reliability (coded units).

Term	Effect	Coef.	SE Coef.	t-Value	p-Value	VIF
Constant		0.89313	0.00312	285.80	0.000
nw (rpm)	−0.01625	−0.00813	0.00312	−2.60	0.032	1.00
nt (rpm)	−0.00125	−0.00062	0.00312	−0.20	0.846	1.00
vf (mmprev)	−0.01625	−0.00813	0.00312	−2.60	0.032	1.00
nw (rpm) × nt (rpm)	0.00125	0.00062	0.00312	0.20	0.846	1.00
nw (rpm) × vf (mmprev)	0.00625	0.00312	0.00312	1.00	0.347	1.00
nt (rpm) × vf (mmprev)	−0.00875	−0.00438	0.00312	−1.40	0.199	1.00
nw (rpm) × nt (rpm) × vf (mmprev)	−0.01125	−0.00562	0.00312	−1.80	0.110	1.00
Regression Equation in Uncoded Units
r = 0.89313 − 0.00813 nw (rpm) − 0.00062 nt (rpm) − 0.00813 vf (mmprev)+ 0.00062 nw (rpm) × nt (rpm) + 0.00312 nw (rpm) × vf (mmprev) − 0.00438 nt (rpm) × vf (mmprev) − 0.00562 nw (rpm) × nt (rpm) × vf (mmprev)
Model Summary
S				R-sq
0.44874				87.31%

and

Table 5. Estimated effects and coefficients for productivity (coded units).

Term	Effect	Coef.	SE Coef.	t-Value	p-Value	VIF
Constant		65,700	300	219.00	0.000
nw (rpm)	5400	2700	300	9.00	0.000	1.00
nt (rpm)	12,600	6300	300	21.00	0.000	1.00
vf (mmprev)	10,200	5100	300	17.00	0.000	1.00
nw (rpm) × nt (rpm)	4200	2100	300	7.00	0.000	1.00
nw (rpm) × vf (mmprev)	1800	900	300	3.00	0.017	1.00
nt (rpm) × vf (mmprev)	9000	4500	300	15.00	0.000	1.00
nw (rpm) × nt (rpm) × vf (mmprev)	−1800	−900	300	−3.00	0.017	1.00
Regression Equation in Uncoded Units.
q = 65,700 + 2700 nw (rpm) + 6300 nt (rpm) + 5100 vf (mmprev) + 2100 nw (rpm) × nt (rpm)+ 900 nw (rpm) × vf (mmprev) + 4500 nt (rpm) × vf (mmprev) − 900 nw (rpm) × nt (rpm) × vf (mmprev)
Model Summary
S				R-sq
0.39170				82.58%

. The third column (in Table 3, Table 4 and Table 5) represents the regression coefficient (a half effect of each factor)

The results of the wood-milling process were observed, and then, dispersions were analysed before and after full-factorial optimisation. The normality of the data distribution allowed for the use of parametric tests. We chose the F-test to verify the significance of variability reduction in terms of the environmental friendliness of production. This method consists of four steps: state the hypotheses, formulate an analysis plan, analyse the sample data and interpret the results. The test criterion is calculated according to:

$$\mathrm{F}=\frac{{\mathsf{\sigma }}_1^2 }{{\mathsf{\sigma }}_2^2}=\frac{{\mathrm{n}}_1\times \left({\mathrm{n}}_2-1\right)\times {\mathrm{S}}_1^2}{{\mathrm{n}}_2\times \left({\mathrm{n}}_1-1\right)\times {\mathrm{S}}_2^2} ,$$

(20)

which has a Fisher–Snedecor distribution of F(n₁ – 1, n₂ – 1).

The null hypothesis is determined by H1: ${\sigma }_1^2$ = ${\sigma }_2^2$.

The alternative hypothesis is then H1a: ${\sigma }_1^2$ ≠ ${\sigma }_2^2$.

If $F>F_{\frac{p}{2}}\left(n_1-1,n_2-1\right)$, we reject the hypothesis H₀₁ (H₁₁ is accepted). In this case, we chose a significance level of p = 0.05. We determined the required characteristics in both groups (swapping the order so that F > 1), and we obtained the following results:

Before full-factorial optimisation:

n₁ = 278; s₁² = 0.6201.

After full-factorial optimisation:

n₂ = 315; s₂² = 0.4073.

After substitution into Equation (10), we obtained:

$$F=\frac{{\sigma }_1^2}{{\sigma }_2^2}=\frac{n_1\times \left(n_2-1\right)\times S_1^2}{n_2\times \left(n_1-1\right)\times S_2^2} =1.5231 \ge F_{0.025}(277;315) = \mathrm{FINV}(0.025;277;314) = 1.2562.$$

(21)

The test criterion exceeded the critical value (1.2562) at 277 degrees of freedom of the first set and 315 degrees of freedom of the second set. Therefore, H1 was rejected. There was a statistically significant difference between the variances; therefore, the factorial optimisation represents a noticeable improvement in the environmental friendliness of the wood-milling process.

After verifying the positive effect of the factorial design on the production variability (and, thus, also increasing its environmental friendliness), we compared the two mean values of the average profit from the set, each with 52 values (weekly averages of the profit from the given product). The test procedure, called the two-sample t-test, was an appropriate method (Student’s t-test). Student’s t-test is a conventional statistical procedure for measuring the significance of a difference of mean [32]. The steps are the same as in the previous F-test. If we can assume б₁² = б ₂² (we also have checked it with the F-test), we choose the test criterion:

$$\mathrm{T}=\frac{{\mathrm{M}}_1-{\mathrm{M}}_1}{\sqrt{{\mathrm{n}}_1{\mathrm{S}}_1^2-{\mathrm{n}}_2{\mathrm{S}}_2^2}}\times \sqrt{\frac{{\mathrm{n}}_1{\mathrm{n}}_2\left({\mathrm{n}}_1+{\mathrm{n}}_2-2\right)}{{\mathrm{n}}_1+{\mathrm{n}}_2}},$$

(22)

which has a Student’s distribution t (n₁ + n₂ − 2).

The null hypothesis is determined as H2: μ₁ = μ₂.

The alternative hypothesis is then H2a: μ₁ > μ₂.

If T | > tp, then we reject the hypothesis H₀ (H₁ is accepted).

In this case, we chose a significance level of p = 0.05. We determined the required characteristics in both groups with the following results:

Before full-factorial optimisation:

n_{1 =} 52; s₁² = 1.26(×10³ EUR²); M₁ = 29.6 (×10³ EUR).

After full-factorial optimisation:

n₂ = 52; s₂² = 1.17(×10³ EUR); M₂ = 35.4 (×10³ EUR).

After substitution into Equation (22), we obtain:

$$\mathrm{T}=\frac{{\mathrm{M}}_1-{\mathrm{M}}_1}{\sqrt{{\mathrm{n}}_1{\mathrm{S}}_1^2-{\mathrm{n}}_2{\mathrm{S}}_2^2}}\times \sqrt{\frac{{\mathrm{n}}_1{\mathrm{n}}_2\left({\mathrm{n}}_1+{\mathrm{n}}_2-2\right)}{{\mathrm{n}}_1+{\mathrm{n}}_2}}=7.306 \ge \mathrm{t}0.05\left(102\right) = \mathrm{TINV}\left(0.05;102\right) = 1.660.$$

(23)

The test criterion exceeded the critical value (1.660) at 102 degrees of freedom. Therefore, H2 was rejected. There was a statistically significant improvement between the mean profit before and after factorial optimisation.

5. Discussion

The coefficient of multiple determination R-Sq(adj) = 79.59%, which indicates that mean profit equation is well suited to the acquired response data. The model can explain the variability to 79.59% of non-negligible interactions. The optimal settings for manufacturing processes using CNC machining are as follows:

• Workpiece speed: 2 rpm;
• Cutting speed: 1500 rpm;
• Axial speed: 0.5 mmpsec.

The sustainability of industrial woodworking lies in the supply of products that meet customers’ needs and requirements, with a low impact on the environment in terms of material and energy consumption. The sustainability of woodworking is part of a whole set of human activities that, taken as a whole, cannot exceed a level of balance with the planet’s capacity. At present, most environmentally friendly production processes are accompanied by lower economic efficiency. Environmentally friendly production also applies to the processes of the woodworking industry, where environmental friendliness is often accompanied by less progressive technology.

In our methodology, we used factorial proposals to minimise product variability (primarily material savings) along with an accompanying set of factors to reduce the average energy consumption (per workpiece). The effect of these savings in variable costs managed to increase the spread between the average revenues and costs, which also led to a higher average profit.

The significance of factorial optimisation was confirmed by a statistical F-test, which showed a reduction in production variability, and a t-test, which showed a decrease in production costs, reduced energy consumption and increased profit margins. Thus, this study shows that it is possible to reduce the impact of wood processing on the environment by striving for optimisation in the setting of operating parameters. If we consider the optimisation of energy consumption concerning the environment in a highly productive CNC centre, it would be appropriate to divide the optimisation into three groups. The first group would be numerically controlled primary and secondary drives. In the second group, energy appliances have pumps for cooling and tool clamping. In the third group, energy appliances have an electrical equivalent of air consumption for CNC machines connected to a common distribution. For the first group of appliances (servo-drives of motion axes), it is possible to reduce energy consumption by modifying the parameters of the cutting process (we focused on this in our article), further lightening the moving materials, and reducing passive resistance (e.g., reducing lubricant viscosity). More energy savings can be made in the area of the recuperation of braking energy back into the system. The energy savings of the second group of appliances focus on the various peripherals of the machines. These include more complex units such as tools or pallet change and fluid management systems. In the third group (components for the distribution and use of compressed air), it is advisable to minimise the pressures and the amount of air used for permanent functions, such as spindle overpressure.

From the above, it is clear that, with modern, highly productive CNC machines, energy saving is a relatively complicated process. Moreover, it is relatively easy to optimise the energy consumption of wood processing production in an older single-purpose machine, which was the subject of our optimisation (Copy Milling Machine COSMA-SKS 7). Our full factorial optimisation (considering the interactions between factors) and tool change has reduced energy consumption per workpiece by 7.2% (see the last row of Table 1). Furthermore, the factorial optimisation reduced production variability by 11.3% (expressed by the relative change in standard deviations delta/s1 = (0.887 − 0.787)/0.887 = 0.113) and increased the average production profit by 16.4%. At present, the company where an optimisation design of the energy consumption of wood processing production was carried out has purchased a highly productive Greda Venus CNC copying machine, Prague [CZ]. For this new device, the optimisation of energy consumption will be performed in the context of savings in the three groups of components mentioned above. Before full factorial optimisation, we created a model of sustainability in the use of resources to define the impact of the optimisation of CNC machines of wood processing production, respecting the sustainability aspect, which is described in Figure 1 and characterised in the article’s conclusion.

To fulfil the theoretical goal, we created a causal description of sustainable resources using the input–output model of a closed system of the wood processing industry, which interacts with the environment. The model serves both as a causal description of the use of renewable and non-renewable sources and as a way of potentially predicting changes in CO₂ concentrations, depending on the amount and ratio of renewable and non-renewable sources. Thus, by solving Equation (8), we predict changes in CO₂ concentrations over time without the need to know the retrospective period and with less sensitivity to changes in environmental factors compared to using time series analysis or trend extrapolation for the regression function. This new method of prediction has a number of advantages over regression extrapolation. The new method does not need to know past data. It just needs to know the current state. Because various influences are included in the initial equation, this prediction method is not burdened by the influence of the third factor. This case occurs when the correlation is apparent and the variable is only an intermediate variable. For example, in Toth et al. [43], the prediction of the number of employees in forestry is based on a regression function with a known time evolution since 1930, i.e., data for the past 89 years. For prediction, these data should be collected for another 44 years, because it is generally recommended that one knows two-thirds of the past for the prediction of one-third of the future [42,44]. Figure 1 (in [43]) shows that, on the one hand, the linear trend is incorrectly described by the quadratic equation; especially in a forecast of six years into the future (2025), the declining regression line extends into the area for negative values of forestry workers. Thus, although the future will behave similarly to the past and thus should have a high correlation coefficient (here it is 0.77), the regression extrapolation still leads to misleading predictions. For this reason, we try to use an alternative mathematical forecasting procedure—in this case, for the prediction of CO₂ in the atmosphere. We have already verified the functionality of this procedure and we plan to publish it in a special issue of Sustainability.

6. Conclusions

Sustainable development is a method for developing a society that reconciles economic and social progress with the preservation of the environment. One of the main goals of sustainable development is to preserve the environment for future generations; from an economic point of view, the main effort is to create a world without consumption, that is, to create an ever-higher quality of products that will enable sustainable prosperity and form a welfare state, moving towards social welfare. Therefore, sustainable development is built on social, economic and environmental pillars that should be considered equally. Currently, sustainable development is, mistakenly, only linked to the environmental pillar. At the same time, both the environmental and social pillars lag behind because of economic growth.

The practical purpose of the present article was to investigate the correlation between the process parameters of (wood) milling and economic variables. We calculated multiple regression equations from common correlations and determined the interactions between individual factors. This progress allowed us to perform an optimisation according to multiple factors and according to significant integrations between factors. The practical objectives of the experiment were two-fold: the first objective was to identify the key milling process parameters that influence the responses of select economic variables. The second objective was to identify the optimal settings for either minimising the energy and unit variable costs or maximising the overall production profit.

These practical objectives were met. First, we established a regression model for productivity, reliability, variable cost and average profit. We found the optimal parametric settings for these variables. We also optimised the setting of parameters for variability and average energy consumption (mean profit). We subjected the obtained values to statistical tests, where we demonstrated the effect of variability and energy savings at a significance level of 0.05.

The limits of the presented solution lie in the statement ‘correlation does not imply causation’, which refers to the inability to confirm a cause and effect relationship between two variables based on an observed correlation. A random factor causes another limitation in determining the parameters of the regression dependence. We eliminated this random factor effect by reducing the production variability by replicating the trials. We also performed randomised experiments to reduce the autocorrelation of the responses. We performed two verifications to reduce the likelihood that the correlation does not imply causality. The first verification was based on the fact that the cause precedes the consequence (the setting of factors precedes the response). The second verification was based on eliminating the influence of the third variable. This second verification was performed using factor screening.

In future research, the authors would like to focus on finding a method for verifying the causality of regression models in the woodworking industry, which were presented here as a factorial design. This research promises to increase the reliability of factorial optimisations in the woodworking area (specifically excluding the influence of the third factor). Moreover, it will bring about a general improvement in the understanding of causal relationships in reducing variability and increasing efficiency in the woodworking industry. The next phase of the methodology is to create a more advanced response surface methodology (RSM) by adding centre points and axial points to the current design.

The whole procedure was carried out on Stusek-DVB Company Ltd., Prague [CZ] as a case study, making it easy to repeat the process. The results of the experiment stimulated the engineering team to extend the applications of full-factorial design to other milling processes for performance improvement and variability reduction.