Saosys toolbox as Matlab implementation in the elastic‐plastic analysis and optimal design of steel frame structures

Abstract The improved mathematical model of steel frame structures’ design is created. The loading is simple, and plastic strains are evaluated. Energy principles of deformable body mechanics and mathematical programming theory are employed. Equilibrium finite elements with interpolation functions of internal forces are used for discretization. The elements are designed using HE, IPE, RHS steel profile assortments and considering dispersion of geometrical characteristics of profile assortment sets. Optimal design of steel structures is realized by using the experimental tool system JWM SAOSYS Toolbox v0.42, which was created by the authors in MATLAB environment. SAOSYS architecture operates with object‐oriented finite elements pseudo‐polymorphously. The possibilities of this system are demonstrated by considering a numerical example of optimal design of industrial building frame with strength and stiffness constraints. The assumption of small displacements is adopted for computations.


Introduction
Structural design evaluating plastic strains allows us to exploit carrying capacity of the structure more effectively and create more economical projects (Majid 1974;Atkočiūnas 1999;Choi and Kim 2002;Kaliszky and Logo 2002;Алявдин 2005). It is worth noticing that the assumptions of the limit equilibrium theory are referred to in many papers (Čyras et al. 2004;Karkauskas 2007). Optimization results under plastic collapse criterion are not decisive in every case because the limit feasibility state of optimal structure can be lost even in the case, when plastic collapse due to excessive inelastic strains and displacements is not achieved (Giambanco et al. 1994;Tin-Loi 2000). Therefore, chronologically, the paper (Kaneko and Maier 1981) is most important for the cases, where elasticplastic structures' optimization with stiffness constraints is analysed.
An optimum criterion of a structure can have energetic or a definite physical meaning, i.e. minimum volume (weight) or minimum cost (Prager 1955;Skaržauskas et al. 2005;Kalanta et al. 2009;Šešok and Belevičius 2008). Engineering development in the field of optimal structural design requires some theoretical and practical knowledge, including the fundamentals of structural mechanics, structural design standards (STR 2.05.08;EN 1993EN -1-1: 2005 and, finally, modern information technologies. Further, not sticking to chronology, we will mention only some papers (concerning structures' optimization) and provide a comprehensive list of literature. It is mostly the work of M. I. Reitman (1976) in Russian, as well as the books (Brandt 1978;Atrek et al. 1984;Lloyd Smith 1990).
It should be noted that the above papers do not take into consideration the relationship between dual theory of mathematical programming and the problems of static and kinematic formulas of rigid deformable body. Namely, the dual theory applied to holonomic plastic deformation process (Koiter 1960) allows the construction of matrices showing the influence of residual internal forces and displacement influence matrices (Atkočiūnas 1994). Finally, the revelation of the mechanical meaning of Kuhn-Tucker optimality constraints (Bazaraa et al. 2004) facilitates numerical realization of optimization problems (Atkočiūnas et al. 2003. An attempt was made to avoid all these imperfections in the current paper. An improved mathematical model of minimum volume design of steel structures with plastic strains was created applying energy principles of structural mechanics and mathematical programming theory. Besides, strength, stiffness and stability requirements to structures discretized by finite elements and subjected to local forces were evaluated more accurately. The extreme internal forces of the elements were also restricted by additionally introduced nonlinear yield conditions. There is also a possibility of precise evaluation of extreme element deflections under stiffness conditions (these conditions are mostly the constraints of node displacements of the discrete model of the structure). These additional means allow us to avoid densification of the finite elements grid, thereby decreasing the size of the optimization problem, as well as saving computer resources (especially, solving time).
An approximate objective function expressing the volume of a structure is used in the developed mathematical model of the optimization problem. Structural elements are designed considering the dispersion of geometrical characteristics in the sets of profiles and based on the principle of admissible fields of geometrical characteristics of the profiles assortments HE, IPE and RHS (Rectangular Hollow Section).
For practical realization of the optimization problem, the authors created a design algorithm of elastic-plastic structures and structural analysis, as well as optimal design system JWM SAOSYS Toolbox v0.42 (Structural Analysis and Optimization System) in MathWorks MATLAB environment (Jankovski and Atkočiūnas 2008) ( Fig. 1). Nonlinear optimization problem considered is nonconvex. The convergence is obtained by an iterative method, i.e. by solving a sequence of nonlinear problems. The system SAOSYS combines finite element method, object-oriented programming (OOP) paradigm, mathematical models (Čyras et al. 2004) based on structural analysis and extreme energy principles of optimization, as well as mathematical programming theory and methods (Bazaraa et al. 2004;Raue et al. 2009), principles of the initial structural data input and parameterization, databases of steel profiles' assortments and the output and interpretation methods of textual and graphical data. The main parts of the optimal design system of structures SAOSYS system's architecture embraces User mode pre-post-processor control system, and Kernel mode, involving system topology and its components. The created system operates with object-oriented LINK11 and BEAM31 finite elements, composing finite elements' library FELIB, which contains equilibrium and displacement formulations of finite elements for solving various problems. Due to limitations of MATLAB 7.0 environment OOP facilities, pseudo-polymorphism is realized in SAOSYS system which maintains a concept of pseudovirtual methods. Pre-processor of the created system employs the command data input method of structural model (similar to ANSYS software) and a possibility of structural parameterization in a variational design case.
The possibilities of SAOSYS system and the proposed technique are illustrated by a numerical example of optimal design of industrial building frame with strength and stiffness constraints. The assumption of small displacements is evaluated for computations.

A mathematical model
An elastic-plastic beam structure of known geometry subjected to specified loading is analyzed. Simple loading is perceived as loading, when all loads are proportional to one common factor: thus, plastic deformation holonomicity is indirectly provided. The principle of the minimum energy of A. Haar and von Kármán Th. (Koiter 1960) is valid for this process. The numerical example in this paper is aimed at finding the project of the structure of minimum volume V, whose optimality criterion (1) is provided with respect to strength and stiffness, as well as stability requirements. The optimality criterion consists of the following items: L is the structural element's length vector; G 0 is the leading vector of the elements' cross-sections of design geometry; A(⋅) is the vector function of cross-section geometry conversion to crosssection areas. Thus, minimization is performed when the whole structure's configuration, physical-mechanical characteristics of the material of the elements, loading F and the vectors of the limiting values of structural nodes' displacements u min , u max and ultimate deflections v min , v max of the elements are given. The vectors v min , v max are used to control stiffness conditions. The basis of this optimization is a system of equations and dependencies, defining a real stress-strain state of an elastic-plastic structure before plastic failure. This system can be created assuming the constraints of the static formulae of the analysed problem and Kuhn-Tucker's optimality conditions of the problem (Atkočiūnas and Merkevičiūtė 2003). The above-mentioned system of dependencies (giving the results of influence matrices [Z] and [Y] of residual internal forces S r and residual displacements u r respectively) is often referred to as a generalized Lagrange problem. Further, the internal forces and displacements at the elastic stage will be denoted as S e and u e , respectively (16, 15), where [β] is the matrix of the influence of structural displacements. Thus, the aim of the computations is to find optimal distribution of geometrical characteristics of the elements' crosssections G 0 , while safe exploitation of the structure with plastic strains is secured.
Thus, a mathematical model for the problem of the minimal volume is as follows: Thus, the mathematical model (1)-(9) for the volume minimization problem of elastic plastic beam structure with variables G 0 and λ consists of nonlinear objective function (1) and constraints-conditions, such as linear inequalities (2, 8); nonlinear complementary slackness condition (3); stiffness constraints (4, 5); additional nonlinear yield conditions (6), evaluating extreme internal forces of elements; nonlinear stiffness constraints (7), evaluating extreme deflections of elements; and constructional constraints (9). Furthermore, the following notation is used in the given mathematical model: Then, the vectors of the total displacements u and internal forces S of the structure are as follows: [ ] It only remains to note that structural configuration matrix [B] and the vector of plastic multipliers λ are included into the direct yield conditions (2) (the vector b is known during the iterative process). The vectors u min , u max are interpreted as maximum vectors of negative and positive values of the restricted nodes' displacements; vectors v min , v max are maximum vectors of negative and positive values of elements' deflections; [A] is the matrix of coefficients of equilibrium equations of the structure; is the stiffness matrix of structural elements; [Φ] is the matrix of structural elements' linear yield conditions at the nodes.

Beam finite elements: yield and strength conditions
Steel structures will be modelled by the equilibrium finite elements with interpolation functions of internal forces (De Veubeke 1963;Gallagher 1975;Kalanta 1995;Wilson 2002). LINK11 and BEAM31 types of finite elements are described in this paper. All elements k = 1, 2, 3, ..., n e of the structure compose a set K of finite elements. Subsets K 11 and K 31 , corresponding to finite element types LINK11 and BEAM31, compose the set K. The subset R r is composed of the elements of the same type, material and cross-section's geometrical characteristics. The set R is composed of the subsets R r of attributes. After describing the sets of the finite elements of structural model, we will discuss yield and strength conditions of every element in detail.

LINK11.
It is an elastic-plastic equilibrium finite element of truss which can only lengthen or shorten (i.e. strains in axial direction are only evaluated) (Fig. 2). The vector of the internal forces of the element is S k .=.{N k } and its yield-strength conditions (i.e. strength and stability requirements) are as follows: The first constraint is the yield condition of the tensioned element section, since the second is the strength-stability condition of the compressed element (STR 2.05.08. 2005). Buckling reduction factor ϕ k of centrally compressed beam is calculated as follows: The yield-strength conditions (20), written in matrixvector form are as follows:  Here, N k is the element's axial force; f y,k is tensile steel strength depending on the yield stress; E k is elasticity modulus of steel; l k is the element's length; A 0,k , I k , i k are the element's cross-section area (design parameter), the inertia moment and inertia radius of cross-section, respectively; [Φ k ], B k (⋅) is the matrix of yield-strength conditions ratios and configuration vectorial function.
BEAM31. It is an equilibrium finite element of the 2D beam under bending and tension or compression (Fig. 3). Two cases of this element are considered: C 1 is an element subjected only to nodal loads (default); C 2 is an element subjected to nodal loads and uniformly distributed load 7 k F ′ . Then, the vector of the element's internal forces S k is as follows: Yield conditions must be satisfied in the nodes of the BEAM31 finite element: , ,0, k p ly k k 1, 2, 3 ; These conditions written in the matrix-vector form are as follows ( ) where M ki denotes bending moments in the element's nodes; ξ.=.0,85 is the ratio; W pl,y,0,k , A k is plastic section modulus of the element's cross-section (design parameter) and cross-section area, respectively. The extreme bending moment M k,ext should be evaluated while designing elements of C 2 case (i.e. when bending moments vary according to the second degree curve) (Fig. 3). It can be implemented approximately by increasing the number of finite elements. The second way of accurate calculation, which will be discussed below, is the direct application of additional nonlinear yield conditions (6). With reference to the formula (19), the true vector of internal forces S.≡.{S k , k.∈.K} can be calculated. A relative position η k of the feasible extreme bending moment is expressed as follows:

Stiffness condition of the BEAM31 finite element
The stiffness conditions (4, 5) in the structural volume minimization problem allow us only to restrict the real (total) displacements u of the nodes of the discrete structural model. In order to check extreme deflections v k,ext ( Fig. 4) of separate elements, we must densify the grid of finite elements, or calculate accurately, i.e. directly apply stiffness conditions (7), corresponding to extreme deflections. We will discuss it in more detail. Applying the conditions (18) we begin to calculate the vector u of global displacements, from which we pick the values of nodes' displacements u k of a single element. 2||3

Fig. 4. Extreme deflections of the BEAM31 element
To create the interpolation function of the element's deflections v k (x), we choose a third-degree polynomial: where [N v,k (x)] is function-matrix of the element's shape; [T k ] is transformation matrix of the displacements of the element's nodes; u k =.{u k1 , u k2 , ..., u k6 , <u k7 > C² } is global displacement vector of the element's nodes; k ′ u is local displacement vector of the element's nodes. For deflection interpolation function v k (x) we apply a stationary condition. Then, while solving the quadratic equation, we try to find two solutions to the locations of extreme deflections x k,ext,i : The extreme deflections of the element are calculated for each location x k,ext,i : Finally, these nonlinear stiffness conditions (which evaluate extreme deflections and are denoted in the mathematical model (1)-(9) as (7)) are expressed as follows:

Assortments: fields of discrete geometrical charac-teristics of the profiles in structural optimization
Steel structures' design is closely connected to discrete sets of profiles' assortments. Analysing the distribution of geometrical characteristics of the profiles IPE, HE, RHS ( Fig. 5, Fig. 6), we can see that the single-valued dependence between cross-section geometrical characteristics A-W pl,y and I-A does not exist. Therefore, the admissible points k G ( Fig. 6) are to be found in discrete fields In the case of the elements of different types in a structure, the geometry vector k G of cross-section takes the form of: For the whole structure, the following notation is correct: The mathematical model (1)-(9) involves only G 0 design parameter (the problem's variable), which is the vector of the leading design geometry. The inertia moments I k and areas A k of cross-sections compose the vector of the driven geometry G 1 . While solving the optimization problem (1)-(9) by the iteration process, the leading geometry is optimized, whereas the lagged driven geometry is only corrected with reference to the yield conditions (2, 6) and admissible field bounds A min (W pl,y

Optimality criterion of the structure
The nonlinear objective function (1) reflects the volume V of the structural elements written using the vector of the leading cross-section geometry G 0 : Functions A k (⋅) relate to different geometrical crosssection characteristics of the elements (G 0,k :.A 0,k .|.W pl,y,0,k ) with the areas A k of these cross-sections. Because of the lack of a single-valued dependency between A(W pl,y ) (the case of the BEAM31 element design) (Fig. 5), we analyzed two appropriate solving methods of this problem: 1..approximate mean curves of profile assortments; 2..isocurve fields, approximating admissible fields of profile assortments A-W pl,y .
The numerical experiments showed that the first method proved to be better because a more optimal solution was reached. The second method 'disbalances' the optimization problem during the iteration process and, therefore, the solution departs from the optimal version. Approximative mean curves. For approximating discrete points' distribution A-W pl,y of the profile assortments HE and IPE we choose a third-degree polynomial ( )  (43) the ratios a i of which are derived by creating the following non-correlation function of the least squares method We apply stationary conditions for the total noncorrelation of a set of discrete points D and a condition for the polynomial, approximating the edge point dependency of a discrete set. Then, four-equation system takes the form of: We solve the given equation system for separate profile assortments at the point of the ratios a i (i = 1, 2, ..., 4). The obtained ratios' values of approximative polynomials (Fig. 7)  A W , we analyzed the fields of the isocurves, approximating the admissible profile fields D (Fig. 8).
The function of isocurves' field is linearly interpolated between two spline functions A min (⋅) and A max (⋅), which bind the admissible field D . This function may be written as follows We solve the optimization problem (1)-(9) by iterations and control the dependence of the k G.≡.{W pl,y,0,k ,.A k } points on the admissible field. During this process, we correct relative ratios of the approximative isocurves (46) selectors: . (47) Then, the optimization problem developed for a new iteration is solved with the new corrected approximative functions , which compose the objective function (1).

Design algorithm
Design algorithm of elastic-plastic structures is realized in the integrated design system, which combines MATLAB and system SAOSYS developed by the authors. Further, we will describe the main parts of the algorithm (Fig. 9).
Structural modelling. Parameterization (in the case of variant design) and discretization of a structure with elastic-plastic equilibrium beam LINK11 and BEAM31 finite elements are performed. The initial data file of the structure (SAOSYS pre-processor Batch and Data file -BDF) is created.
Preparing design environment. The routine P1 of SAOSYS pre-processor reads and translates BDF, as well as creating the database DB of structural finite element model. The routine P2, which controls profile assortments, reads and prepares the assortment of steel profiles HE, IPE, RHS from SRT database. The routine P3 loads the finite elements' library FELIB of the SAOSYS system. Finally, the routine P4 creates and initializes the assemblage FE of the finite elements, which compose the structure (i.e. calculates the length of the elements l k and direction cosine vectors n k , as well as preparing displacement compatibility vectors and performing other operations).
Preliminary calculations. The routine P5 directly creates the ratios' matrix [A] of the structure from separate finite element matrices [A k ] of assemblage FE and external loading matrix [F] (in the case of single loading -vector F). The boundary conditions of a discrete model are evaluated. The routine P6 activates the yield conditions of the elements of FE assemblage. The edge values' vectors G 0,min , G 0,max , G 1,min , G 1,max of the leading and driven geometry with reference to SRT database of profile assortments are created. The total length vector L of the sets of structural elements and element length vector L max composed of the longest elements in groups (R r , r.∈.R) are created. The vectors of the nodes' limit displacements u min , u max and the element's ultimate deflections v min , v max are also created.
Solving the optimization problem. To solve the optimization problem (1)-(9), we use an iterative approximation and begin with the highest geometrical values of the vectors G 0 .=.G 0,max , G 1 .=.G 1,max .
•_Step_1: the routine P7, with reference to the created geometrical matrix [G] (41) of W pl,y,0,k and A k values, performs an interpolation operation I y,k .=.I y (⋅) of the cross-sections' inertia moments of the BEAM31 finite elements. This operation is described in more detail in Section 7.1. The created vector I y of interpolated inertia moments is used for constant flexibility matrices [D k ] for BEAM31 elements.
•_Step_2: matrices and vectors B k (23, 27) finite elements' flexibility [D k ] (Gallagher 1975;Kalanta 1995) and yield conditions [Φ k ] are created: •_Step_3: with reference to formulas (15, 16), the displacement vector u e and internal forces' vector S e for the elastic solution of the stress-strain state of the structure are found.
•_Step_4: the routine P8 solves one iteration of nonlinear mathematical programming optimization problem (1)-(9). If the procedure of solving is successful (i.e. optimal solution is found), we have a new leading geometry vector * 0 G of cross-sections and a new vector.λ of plastic multipliers. (1)-(9) Evaluate: u, S

Correct G1
Evaluate Volume V

Is Convergence
Achieved?
Optimal Solution: Step 6 Step 7 Step 1 Step 5 Step 4 Step 3 Step 2 dLINK1 dBEAM3 If the solution failed (i.e. the admissible point and optimal solution are not found), the leading geometry vector G 0 is increased: o t h e r w i s e , is the leading geometry vector of previous successful iteration; ε is the relative threshold of recurrent G 0,k increase (10 -3 %); ξ is the partial ratio of G 0,k direct increase. The routine P9 corrects the driven geometry vector of cross-sections G 1 ; then, we return to Step_1.
•_Step_5: with reference to formulas (18, 19) the displacement vector u and internal forces' vector S of the real (total) stress-strain state are calculated.
•_Step_6: the routine P9 performs a correction procedure of the driven geometry vector G 1 . The procedure is described in detail in Section 7.2.
•_Step_7: with reference to geometry matrix [G * ] cross-section areas 0,k A * and k A * values, the structure's volume V is calculated. This iterative process is performed until the above convergence conditions of the problem are satisfied: where V′ is structural volume of the previous iteration; 0 , G V ε ε denote convergence tolerance criteria (0,1 %) of the leading geometry of cross-sections and structural volume, respectively. The results of problem solution verification. Based on the optimal solution * * 0 1 , G G , u, S, V, the postprocessor of the system SAOSYS calculates the distribution of strength reserves in the elements' length and cre-ates the control diagram EYCPlot of the elements' strength reserve. In addition, the internal forces and structural strain intensity diagrams can be created, and the numerical results of node displacements, as well as extreme deflections and internal forces of the elements can be derived.

Interpolation of the inertia moments of cross-sections
Let us assume that the admissible point of cross-section k G belongs to the discrete field which can be written in the matrix form as follows and, finally, linear interpolation ratios can be derived from this formula 1 1 2 (56)

Correction of cross-section geometry G 1
In structural design, we introduce a concept of the set of the elements' subsets R. We optimize cross-sectional geometrical characteristics G 0,r , G 1,r of separate subsets of elements r.∈.R. These characteristics compose a couple of vectors G 0 and G 1 . Since we operate with the subsets of elements, the elements of the subset r are obtained by the intersection of element sets K and R r (K.∩.R r ). The driven geometry vector G 1 of cross-sections is treated as the limit geometry vector, which satisfies yield conditions (2, 6) of the elements and the limits of admissible discrete fields , r k K R r R ∈ ∩ ∈ .
LINK11. We will deal with estimation of the inertia moments I lim,r of LINK11 elements' limit cross-sections (57), which satisfy the yield conditions (20) and discrete limits of assortments. Let us note that A 0,r .≡.G 0,r . For every compressed element (N k .<.0) of the set R r of the elements' subsets we calculate limit buckling ratios as follows: While calculating the limit buckling and slenderness ratios of this system of conditions: It is additionally checked and corrected. Then, in the cases of the elements under compression or tension, the limit conditions of discrete admissible fields of profile assortments must be satisfied ( ) ( ) 0, 0, min,r r lim,r max,r r Since ( ) ϕ ⋅ is a gradually decreasing piecewise function, the binary search algorithm realizes a numerical solution of the system of conditions (59) (Fig. 11). This algorithm can be implemented by using MATLAB-SAOSYS routine x = BinFunArgValSearch(hFun, y, vInt, tol), where x is the function argument value found (λ lim,r,k ); hFun is the function handle; y is the function output (ϕ lim,r,k ); vInt is the search interval vector {0;.λ c,u }; tol is the search tolerance (convergence tolerance criterion). where vp is the index vector of discrete points found in the bar; vD is the vector ( D A HE∨IPE∨RHS ) of discrete values arranged in the increasing order; x is the real value (A 0,r ); b is the width of the search bar. According to the vector vp, we get discrete limits and return the non-admissible points r G.=.{A 0,r ,.I lim,r } (Fig.12, points 1 and 3) to the admissible field (points 4 and 5).
BEAM31. We will deal with estimation of limit cross-section areas A lim,r (57) of the BEAM31 elements which satisfy the yield conditions (25, 32) and discrete limits of assortments. Let us note that W pl,y,0,r .≡.G 0,r . Similarly, for every element k.∈.K 31 .∩.R r , we calculate the limit cross-section areas A lim,r,k , which compose the vector A lim,r of the limit cross-section areas: where , , , , ,0, Finally, we get the limit cross-section area of the set of the elements' subsets R r as follows: It is also additionally checked and corrected. Then, the following limit conditions of discrete admissible fields of assortments must be satisfied SAOSYS system architecture is based on User mode pre-post-processor control system and private active Kernel mode block, which embraces system topology and its components (Fig. 13). Processors give SAOSYS control commands, functions and scripts to user for modelling and solving the problem and for interpreting textual and graphical data. Further, the development of the experimental system will be aimed at realizing the control of all processor commands. The system's kernel is composed of the library FELIB of object-oriented finite elements; the database SRT of steel profile assortments; the processor modules PRCMDS for problem solving; the database DB of the problem and the structure; the assemblage FE of finite element objects composing the structure; the collection UTILS of additional neutral functions; and the system's functional core KERNEL composed of two subkernels DSK and ESK, working with displacement and equilibrium finite elements' formulations, respectively, for particular problem solving.
Pre-processing module. This system module is intended for the preparation of SAOSYS environment and structural modelling by finite element method, including the creation of load cases. The system pre-processor functions according to the deterministic finite state machine (DFSM) principle for parsing the formatted strings with comma-separated values (CSV) (Aho et al. 1986;Hopcroft et al. 2001;Bucknall 2001). There are two command groups of specifying and action. Specifying commands define the operations, which further are performed by action commands. Command arguments are characterized by the required (noted as arg) and default (noted as arg) values (Table 2).
For the discussed structural discrete model, we perform problem modeling by SAOSYS pre-processor commands, i.e. create Batch and Data File (BDF). The activated pre-processor (using BDF) creates SAOSYS structural model database DB, prepares steel profile assortment database SRT, finite elements' library FELIB, and, according to the data in DB, object-oriented finite elements' model assemblage FE of the structure and initializes the finite elements. Further, pre-processor leaves the work to the selected processor modules. Defines the main problem title.

/EFORM, ef
Chooses a formulation of system finite elements: displacement or equilibrium finite elements.

NDOF, ndof
Chooses the problem type: truss or plane problem.

MP, id, optpars
Defines a list of physical properties of the material.

R, id, optpars
Defines a group of the element's attributes. Processing modules. In the present system version (v0.42), the following PRCMDS processing modules are realized for problem solving: StatAn is the static structure's analysis; EPSOptim is the optimal elastic-plastic steel structure's design (analyzed in this paper); TrussDPD is the direct probability design of optimal steel trusses (Jankovski and Atkočiūnas 2008).
Every processing module has the main subroutine, which prepares a supportive environment of algorithm and performs problem solving. When the solving procedure is completed, the main subroutine copies all the data from the local memory stack (the routine internal workspace) and pastes it into the global MATLAB base workspace memory for further post-processing interpretation.
Post-processing module. This module is intended for output and interpretation of textual and graphical results. The following graphic functions are noteworthy: NPlot() displays nodes; EPlot() is used for deformed and non-deformed structural schemas of displays; SPlot() creates diagrams of internal forces; EYCPlot() creates reserve diagrams of the elements' strength; UPlot() creates the deformed schema of the structure, according to intensity values of displacement u X , u Y and u Σ .

Model database.
All information about the problem and finite elements model of the structure is placed in the SAOSYS database DB of structural model (Fig. 14). It allows us to create a usable structure of the initial data, serving as a basis for finite element model creation. The database consists of the following tables: materials of structural elements MAT; attributes of structural element groups REAL; the discrete structure model nodes NODE; finite elements of the structure ELEM; groups of external loads LOAD; load cases LVAR. The problem is additionally defined by the finite elements' formulation EFORM, degree of freedom of discrete model node NDOF, problem title TITLE and other parameters. Database of steel assortments. The design of steel structures is dependent on databases of profile assortments. Relational Database Management System (RDBMS) and its internal integration to SAOSYS system are refused. Instead, for simplicity, the plain text files, binary files and MATLAB structural principles are used. The system operates with HE, IPE and RHS steel profile assortments. We will discuss the creation of IPE assortment prior to its use in design.

REAL
Steel profile characteristics are written in a not formatted plain text file IPE.sec (Fig. 15). The assortment file has two directives: #srt is the assortment name; #par denotes the labels of column data. Profile names and geometric characteristics are given in the rows below.  Fig. 15. A fragment of profile assortment file IPE.sec SAOSYS function Sortiment() compiles text file IPE.sec and creates binary assortment file IPE.seb, which consists of a header and data segment. Further, SAOSYS works only with binary files. Therefore, assortment data reading into MATLAB workspace memory is fast. Function Sortiment() works in the reading mode. It reads the selected binary assortment and returns data structure IPE, which is defined by the assortment type TYPE; name srtName; names of profile characteristics' parName; and the number of profiles and arrays of characteristics' values. All profile assortments are integrated into SRT database, i.e. assortment table (Fig. 16), which is later used by SAOSYS design algorithms.

Object-oriented finite elements and pseudo-polymorphism
The system of structural modelling, analysis and design SAOSYS operates with object-oriented finite elements. The main choice of object-oriented programming (OOP) was determined by the concepts of encapsulation and polymorphism (Gamm et al. 1995;Riel 1996;Eckel 2000). The library FELIB of the system's finite elements consists of the following classes of finite elements: elastic elements LINK1, BEAM3 and elastic-plastic elements LINK11, BEAM31. Finite element class is a collection of properties (variables) and methods (functions) working with these variables. New finite elements can be integrated into SAOSYS for new problem solving.
An element class constructor creates the finite element object in MATLAB environment memory. Object characteristics and pointers to methods (in nested function forms) are placed inside the object, i.e. in its data field. Standard MATLAB 7.0 object control methods subsasgn() and subsref() are overloaded, and this allowed us to create a compact syntax form, similar to that used in C++ programming language: obj.member = expr The MATLAB functions get() and set() are generally intended for object manipulation (Krysl 2005;Register 2007). Having disposed of these inconvenient functions, we achieved a marked improvement in the program's text clarity and application simplicity. Direct assignment of expression expr result to object member obj.member is performed by (66); the assignment of the results to arrays and cell array members is performed by (67, 68); direct assignment of object member value is performed by (69); assignments of array and cell array member values are performed by (70,71). Syntactic collision of calls to object methods (72) and assignment of object submember (70) is solved by using the lists of method names placed in the classes. These lists separate object members-data from object members-methods.
A polymorphism concept is used to perform flexible operations with finite elements. MATLAB 7.0 system provides limited OOP facilities, not supporting polymorphism. Therefore, pseudo-polymorphism is implemented in SAOSYS system. In this concept, the basic class of finite element represents the cell array FE of structural elements' assemblage. The cell array FE maintains a concept of pseudo-virtual methods (Fig. 17) in the derived classes of elements. The base class of finite elements FE is more oriented to beam elements and has the following properties: dimension degree SPACE; type identificator TYPE; element number ID; the number of the members SNUM of the vector of the element's internal forces S k ; the number of the element's displacement vector u k members UNUM; the number of the yield conditions PHINUM of the element nodes; the element's material identificator MID; the element group identificator RID; the element length L, direction cosines vector vDir; etc. The base class of finite element implements such pseudovirtual methods: element initialization Init(); initiatorselector of the element's yield conditions InitYC(); the getting methods of the transformation matrix [T k ] GetTk(), the matrix of coefficients of the structure's equilibrium equations [A k ] GetAk(), flexibility matrix [D k ] GetDk(), yield conditions [Φ k ] and B k GetYC(); the method of deformed element and evaluation of displacement intensity distribution EvalDefU(); the printing method of the element's internal data EPrint(); the method of the deformed and undeformed element plotting EPlot(); methods of the diagram output of the element's internal forces SPlot() and strength reserve EYCPlot(). The classes eLink11 and eBeam31 of the derived elements are complemented with additional characteristics and overloaded methods, including private helper routines.
First, we choose the formulation of the finite element method (equilibrium or displacement finite elements). Then, the SAOSYS system calls finite element class constructors (eBeam31(), eLink11()) in turn and prepares the library FELIB (the table of objects of the finite element types) (Fig. 18). In MATLAB environment, the array FELIB represents object samples.
While parsing BDF, the SAOSYS pre-processor places finite elements in the element assemblage array FE. Here we can indicate three new element placement steps: 1. selection -queried finite element object selec-tion from FELIB; 2. copying -queried element object copying into FE; 3. initialization -sample element initialization method Init() call, which initializes finite element properties and method pointers.  The array FE, validated by pseudo-virtual methods, prompts to the realization of objects independent of an object type. For example, to derive the matrix of coefficients of the structure's equilibrium equations of k-th finite element, we need the prompt mAk.=.FE{k}.GetAk(). MATLAB 7.0 environment does not allow us to prompt directly to the memory addresses of the data fields (because it does not maintain a pointer data type). Therefore, the values of the variables are updated by using reinitialization principle, which lowers the efficiency of MATLAB environment, especially, in the case of objects. Therefore, contrary to other high-level programming languages (C++, Delphi, Java), it is necessary to update the object by reassigning FE{k}.=.FE{k}.Init(id).

Graphical representation of displacement intensity of the deformed element
The optimization problem (1)-(9) is solved. Then, according to (18), the real displacement vector u of the structural nodes is computed. SAOSYS post-processor routine UPlot() allows us to create the deformed schema, which shows the distribution of displacements' (u X , u Y , u Σ ) intensities in the structure.
Let us discuss a procedure of the finite element BEAM31 presentation (Fig. 19). The vector u k of the finite element nodes' displacement from the vector u is selected. Invoking the function-matrix [H k (x)] for producing the element displacements' interpolation form, we created the element displacements' interpolation vector function u k (x).=.{u k,x (x), u k,y (x)} in the system of the local coordinates (LCS) x1y.
We divide an element into N equal sections (Fig.  19 a). Then, the coordinate vectors P i .≡.{P i,X , P i,Y } in the global system of the coordinates (GCS) are calculated as follows: where [T k ] is the transformation matrix of node displacements of the element GCS→LCS; s is the scale ratio; [ ] k T ′ is the transformation matrix of the coordinates (2×2); J 1 is the vector of the element's first node coordinates in GCS. Let the element be described by the thickness t (Fig.  19 c, d). Then, this condition (the sum of geometric conditions of the vectors (75)) is valid for calculating the vectors of the coordinates L i and R i of the element's multiline nodes , ; We solve the system of equations (75) in respect of the unknown ratios γ and ξ. Finally, the coordinates of multiline nodes are calculated by the formulas given below: The sections of the deformed element R i-1 -R i -L i -L i-1 are covered with quadrangles (Wright and Lipchak 2004), the nodes' colour vectors of which c j .≡.{R j , G j , B j }= c j (u j,X∨Y∨Σ ) are displacement intensity functions of nodes P i in respect of all displacements of structural elements.

A numerical example
Design structure. Industrial building frame subjected to simple loading is designed at the elastic-plastic stage (Fig. 21) Batch and Data File. The preparation of the design system, structural parameterization and sample fragment of modelling commands in pre-processor mode are contained in the file Frame.m of the form given in (Fig. 20). It is a plain text MATLAB script file, beginning with the HEADER consisting of SAOSYS environment preparation macros. Further, ordinary MATLAB computations (calculations of node coordinates, etc.) can be performed. The declared user-named variables (b1, q0) are the structural parameters. They are directly applied to SAOSYS preprocessor's command fragment, which is bracketed by "%{.….%}" symbols (the MATLAB block comment).
The following parameters are defined in the segment of SAOSYS pre-processor commands: the problem title /TITLE; finite element formulation /EFORM (EF_E denotes equilibrium finite elements); the degree of freedom of the node NDOF (ND_PLANE is the 2D problem with three freedom degrees of the node); element materials MP (here, elasticity modulus E and the tensile steel strength depending on the yield stress Ry should be defined); finite elements' sets R (here, sec is the cross-section type, sid is the default cross-section number in assortment and sel is the cross-section selection flag). Further, the nodes' coordinates N of the structure's model are defined. Finite element definition (E command) is performed only after specifying the pointers to the type of the element TYPE, material MAT and the element's attribute group REAL. Then, boundary conditions of nodes (supports) D, and the element's node releases (hinges) ER should be defined. Finally, load cases LOAD are created, where F and SFE commands are assigned the values of concentrated and distributed loads, respectively. At the end of the file, a prompt to SAOSYS processor module routine is written: EPSOptim() is the elastic-plastic structure's optimization with stiffness constraints. Then, the BDF can be executed.
The results obtained. Structural design was performed by using an iterative procedure. In general, eleven iterations were made (Fig. 22). As a result, optimal theoretical cross-sections were found. The profiles closest to them are presented in the table (Table 3). While using the alternative design, it is advisable to choose the frame's second-floor column E{4, 5} cross-sections from the IPE type profiles, while, the first-floor column E{1, 2, 3} cross-sections are chosen from HE type profiles. The volume of the designed structure is V.=.6,897⋅10 -2 m 3 . In the future, we intend to realize discrete optimization of structures, which is required for correcting the final findings of the optimal discrete solution.