Open access peer-reviewed chapter

Constraint Handling Optimal PI Controller Design for Integrating Processes: Optimization-Based Approach for Analytical Design

By Rodrigue Tchamna and Moonyong Lee

Submitted: August 18th 2017Reviewed: January 23rd 2018Published: September 12th 2018

DOI: 10.5772/intechopen.74301

Downloaded: 201

Abstract

This chapter introduces the closed-form analytical design of proportional-integral (PI) controller parameters for the optimal control subjected to operational constraints. The main idea of the design is not only to minimize the control performance index but also to cope with the constraints in the process variable, controller output, and its rate of change. The proposed optimization-based approach is examined to regulatory and servo control of integrating processes with three typical operation constraints. To derive an analytical design formula, the constrained optimal control problem in the time domain was transformed to an unconstrained optimization in a new parameter space associated with closed-loop dynamics. By taking the advantage of the proposed analytical approach, the optimal PI parameters can be found quickly based on the graphical analysis without complex numerical optimization. The resulting optimal PI controller guarantees the globally optimal closed-loop response and handles the operational constraints precisely.

Keywords

  • constrained optimal control
  • industrial PI controller
  • analytical design
  • constraint handling
  • integrating process
  • optimal servo and regulatory control

1. Introduction

Many units used in the chemical process industry, such as heating boilers, batch chemical reactors, liquid storage tanks, or liquid level systems, are integrating processes in which the dynamic response is very slow with a large dominant time constant. In modern control, the integrating process also appears in many applications including space telescope control systems, lightweight robotic arms, and pilot crane control systems. Constraints are inherent in any industrial control systems, either implicitly or explicitly. They are generally associated with both the process variable and controller output. Indeed, typical operational constraints usually include the actuator magnitude and its rate saturation, process/output variable, and internal state variables. The objective of constrained optimal control is to minimize the control cost subjected to constraints on state variables and/or output variables. The importance of taking constraints into account during the design stage of the controller is no more questioned. In fact, a well-designed optimal control would fail in a real-life situation if the constraints are not taken into account while designing the controller. However, optimal control of a process with multiple constraints is still challenging even for a process with simple dynamics. In a popular approach using Pontryagin’s principle or the Hamilton-Jacobi-Bellman equation for a classical optimal control framework [1, 2], the optimal controller parameters are obtained via numerical solution of the nonlinear constrained optimization. However, the existing numerical methods neither guarantee a global optimal solution nor provide useful insights and physical interpretations of the complex relationships existing between the process parameters and control performance. To address this issue, the analytical solutions of optimal proportional-integral (PI) controller under constraints were previously proposed using the optimization-based approach for integrating systems [37] and extended to first-order systems [811]. This chapter introduces the optimization-based approach for the analytical design of optimal PI controller parameters for integrating processes without violating the operational constraints under a unified framework.

2. Formulation of constrained optimal PI control problem

Figure 1 presents the schematic diagram of an integrating process considered in this chapter. It is a type-C PI controller, also called I-P controller, which is a modified type of PID controller where the set point is removed from the proportional term in order to avoid the initial quick on the manipulated variable for a step change in the set point.

Figure 1.

Block diagram of the feedback control of integrating process.

The major resulting transfer functions of this closed-loop system are expressed as

Ys=1τcτIs2+τIs+1Ysps+KpτcτIsτcτIs2+τIs+1DsE1
Us=1KpsτcτIs2+τIs+1YspsτIs+1τcτIs2+τIs+1DsE2

where

τc=1KpKc;Kp=KτE3

The closed-loop damping ratio of the above system becomes

ζ=12τIτcE4

The goal of a constrained optimal problem is to minimize the weighted sum of the process variable error, e(t), and the rate of change in the manipulated variable, u′(t), for a given step change, ΔD/s, in disturbance (i.e., optimal regulatory control) or that, ΔYsp/s, in set point (i.e., optimal servo control) subjected to the following three typical operational constraints: the maximum allowable limit in (1) the controlled variable, ymax, (2) the rate of change in the manipulated variable, umax, and (3) the manipulated variable, umax.

Consequently, the constrained optimal control problem is formulated as

minΦ=0ωyytyspt2+ωuut2dtE5a

subject to

ytymaxE5b
utumaxE5c
utumaxE5d

Through some mathematical operations, the above optimal control problem formulated in the time domain can be transformed to the form with the two new design parameters ζand τcas expressed in Table 1. As shown in Table 1, the three constraints are also expressed as only a function of ζand τc. Then, a simple graphical examination of the contour of the objective function and the constraints in ζτcspace allows to find the location of global optimal solution without complex numerical optimization process.

Regulatory controlServo control
Objective functionminΦrζτc=ατc3ζ2+βτc14ζ2+1minΦsζτc=ατc4ζ2+1+βτc3ζ4
Constraintsτcgrζγghrζτcγhfrζτcγfgsζγghsζτc2γhfsζτcγf
Parametersα=2ωyKpΔD2;β=ωu2ΔD2;γg=ymaxKpΔD;γh=umaxΔD;γf=umaxΔDα=ωyΔYsp22;β=ωu32ΔYspKp2;γg=ymaxΔYsp;γh=KpumaxΔYsp;γf=KpτcumaxΔYsp
Functions g,h,fgrζ=21+x2exptan1xxfor0<ζ<1=2exp1forζ=1=21x2exptanh1xxforζ>1gsζ=1+expπxfor0<ζ<1=1forζ1
where x=1ζ2ζfor0<ζ1=ζ21ζforζ>1
hrζτc=12τcζexp1xtan14ζ214ζ23xforζ<12=u0ΔD=1τcforζ12hsζ=14ζ2forζ>0
frζτc=1+exp1xtan12ζ22ζ21xπxfor0<ζ<12=1+exp1xtan12ζ22ζ21xfor12ζ<1=1+exp2forζ=1=1+exp1xtanh12ζ22ζ21xforζ>1fsζτc=12ζexp1xtan1xfor0<ζ<1=12exp1forζ=1=12ζexp1xtanh1xforζ>1

Table 1.

Objective function and constraints of the optimal control problem in ζτcspace.

3. PI controller design

3.1. Optimal regulatory control

Applying the Lagrangian multiplier [12], it converts the constrained optimization problem in Table 1 to an equivalent unconstrained problem. In regulatory control, the constrained problem can be converted as

minLζτcϖσ=Φrζτc+ϖ1γhhrζτcσ12+ϖ2γggrζτcσ22+ϖ3γffrζτcσ32E6

where ϖiand σiare the Lagrange multiplier and the slack variable, respectively.

The necessary conditions for an optimal solution are then

Lτc=Φrτc+ϖ1hrζτcτc+ϖ2grζ+ϖ3frζτcτc=0E7a
Lζ=Φrζϖ1hrζτcζϖ2τcgrζζϖ3frζτcζ=0E7b
Lϖ1=γhhrζτcσ12=0;Lϖ2=γgτcgrζσ22=0;Lϖ3=γffrζτcσ32=0E7c
Lσ1=2ϖ1σ1=0;Lσ2=2ϖ2σ2=0;Lσ3=2ϖ3σ3=0E7d

The simultaneous solutions of Eqs. (7a)(7d) for possible combinations of σi=0, σi0, ϖi=0, and ϖi0are associated with the corresponding optimal cases. Note that instead of introducing the slack variables, Karush-Kuhn-Tucker conditions [13] can also be utilized for solving the constrained optimization problem, which finds the same optimal PI parameters by the Lagrangian multiplier method.

Figure 2 presents seven possible cases for the location of global optima: the global optimum can be found inside the feasible region (case A), or on the boundary of one constraint (cases B, C, and E), or on the intersection point of two constraints (cases D, F, and G).

Figure 2.

Contours, constraints, and possible locations of the global optimum in regulatory control case.

The global optima of the seven cases can be evaluated by inspecting their geometrical characteristics in ζτcspace as well as the corresponding conditions of the Lagrange multipliers and slack variables as follows:

Case A ϖ1=ϖ2=ϖ3=0: The extreme point, ζτc, which is located inside the feasible region, is therefore the global optimum. Solving Eqs. (7a) and (7b) simultaneously, the global optimum can be determined in explicit form as

ζ=12E8a
τc=βα14E8b

Case B σ1=ϖ2=ϖ3=0: The global optimum, symbolized as ζhτch, is positioned on the constraint, γh=hrζτc. ζhand τchcan be obtained by substituting σ1=ϖ2=ϖ3=0into the equation of necessary conditions, thus solving the following system of equations:

γhhrζτc=0E9a
ΦrζΦrτchrτc1hrζ=0E9b

Case C ϖ1=σ2=ϖ3=0: The global optimum, ζgτcg, is located on the constraint, γg=τcgrζ, and obtained by solving the following system of equations:

γggrζτc=0E10a
ΦrζgrτcΦrτcgrζ=0E10b

Case D σ1=σ2=ϖ3=0: The global optimum represented by ζghτcghis located on the intersection point by γh=hrζτcand γg=τcgrζ, thus can be calculated by solving

γgτcgrζ=0E11a
γhhrζτc=0E11b

Case E ϖ1=ϖ2=σ3=0: The global optimum, ζfτcf, is located on the constraint, γf=frζτc, and can be found by solving

γffrζτc=0E12a
frτcΦrζfrζΦrτc=0E12b

Case F ϖ1=σ2=σ3=0: The global optimum, ζgfτcgf, which is on the intersection point created by γf=frζτcand γg=τcgrζ, is calculated by solving

γffrζτc=0E13a
γgτcgrζ=0E13b

Case G σ1=ϖ2=σ3=0: The global optimum, ζfhτcfh, is located on the intersection point of the constraints γf=frζτcand γh=hrζτc, and calculated by solving

γffrζτc=0E14a
γhhrζτc=0E14b

After the global optimum is determined in ζτcspace, the optimal PI parameters corresponding to each case can then be calculated from Eqs. (3) and (4) as

Kcopt=1Kpτcopt;τIopt=4ζopt2τcoptE15

One of main advantages of the optimization-based graphical approach is that the conditions for the seven possible cases can be directly evaluated based on a meticulous analysis of the graphical shape of the constraints and contours in ζτcspace. The concept of the relative locations between the extreme point and its projections to the constraints is used mainly to develop the conditions to discriminate each case associated with the corresponding global optimum. Figure 3 shows an example of the projection of the extreme point and its notation rule used in this graphical analysis for the optimal PI design. The notation, ζf, represents the abscissa of the projection of the extreme point on the constraint curve by f. Similarly, τchand τcgindicate the ordinate of the projection of the extreme point on the constraint curves by hand g, respectively. If τcis such that τchτcτcg, then the global optimum is above the constraint, h, and below the constraint, g. Referring to Figure 2, this corresponds to cases A, E, or possibly G. In this case, it is apparent from Figure 2 that if ζf<ζ, it belongs to case A (i.e., the extreme point is the global optimum), otherwise it belongs to either case E or G. Cases E and G can be distinguished simply by comparing τcfand τcfh, where τcfhis the ordinate of the intersection of the constraints by fand h. As seen in Figure 2, if τcf>τcfh, the global optimum case belongs to case E, otherwise case G.

Figure 3.

Projection of the extreme point on the constraints in ζ τ c space.

Using similar reasoning, the conditions to discriminate each of seven cases associated with the global optimum can be established according to the relative locations between the extreme point and its projections to the constraints. Table 2 lists the results for the conditions and characteristics of the global optima.

CaseConstraint specificationConditionGlobal optimumLocation of global optimumCalculation of global optimum
AMild ymax
Mild umax
Mild umax
ζfζτchτcγggrζζτcIn the interior of the feasible regionζ=12τc=βα14
BMild ymax
Tight umax
Mild umax
ζfζhζghτc<τchζhτchOn γh=hrζτcγhhrζτc=0ΦrζΦrτchrτc1hrζ=0
CTight ymax
Mild umax
Mild umax
ζfζgζghτcfτcgfτcγggrζζgτcgOn τc=γggrζγggrζτc=0ΦrζgrτcΦrτcgrζ=0
DTight ymax
Tight umax
Mild umax
ζhζghτc<τchor ζg>ζghτcγggrζζghτcghOn the vertex by τc=γggrζandγh=hrζτcγgτcgrζ=0γhhrζτc=0
EMild ymax
Mild umax
Tight umax
ζ<ζfτcfhτcfor ζgh>ζgτcgfτcfτcfhζfτcfOn γf=frζτcγffrζτc=0frτcΦrζfrζΦrτc=0
FTight ymax
Mild umax
Tight umax
ζf>ζg>ζghτcf>τcgfζgfτcgfOn the vertex by γf=frζτcandτc=γggrζγffrζτc=0γgτcgrζ=0
GMild ymax
Tight umax
Tight umax
ζ<ζfτcfτcfhor ζghζhζfτc<τchζfhτcfhOn the vertex by γf=frζτcandγh=hrζτcγffζτc=0γhhζτc=0

Table 2.

Global optima of the constrained optimal regulatory control problem.

3.2. Optimal servo control

The constrained optimization problem in Table 1 can be converted into an equivalent unconstrained problem by applying the Lagrangian multiplier [12] as follows:

minLτcζϖσ=ατc4ζ2+1+β4τc2ζ2+τ2τ2τc3ζ4+ϖ1τc2γhhsζσ12+ϖ2γggsζσ22+ϖ3γffsζτcσ32E16

where ϖiand σiare the Lagrange multiplier and the slack variable, respectively.

Applying the same way used in the regulatory control case, the seven optimal cases can be found by solving the necessary conditions of the above unconstrained problem for the corresponding combination of slack variable and Lagrange multiplier. Figure 4 illustrates the seven possible locations of the global optimum.

Figure 4.

Contours, constraints, and possible locations of the global optimum in servo control case.

After obtaining the global optimum for a particular case, the optimal parameters of the PI controller can be calculated using Eq. (15), i.e., Kc=1/Kpτcopt;τI=4ζopt2τcopt. Table 3 summarizes the conditions that lead to each global optimal location.

CaseConstraint specificationConditionGlobal optimumLocation of global optimumCalculation of global optimum
AMild ymax
Mild umax
Mild umax
ζminζτcmaxτchτcf
where ζminis the minimum allowable damping factor by solving gζmin=γg
ζτcIn the interior of the feasible regionζ=12τc=4βα14
BMild ymax
Tight umax
Mild umax
ζh>maxζhfζminτc<maxτchτcfζhτchOn τc2=γh1hsζζh=4βγh2α+1412τch=γh1hsζh12
CTight ymax
Mild umax
Mild umax
ζmin>ζτcg>maxτcghτcgfζgτcgOn gsζ=γggsζg=γgτcg=1ζg3βα4ζg2+114
DTight ymax
Tight umax
Mild umax
τcmaxτchτcfζmin>ζτcgh>maxτcgτcgfor τc<maxτchτcfζgh>maxζhfζhζghτcghOn the vertex by gsζ=γgandτc2=γh1hsζgsζgh=γgτcgh=γh1hsζgh12
EMild ymax
Mild umax
Tight umax
τc<maxτchτcfζmin<ζf<ζhfζfτcfOn fsζτc=γfγf=fsζfτcffsζτcτcΦsζfsζτcζΦsτc=0
FTight ymax
Mild umax
Tight umax
τcmaxτchτcfζmin>ζτcgf>maxτcgτcghor τc<maxτchτcfζf<ζmin<ζhfζgfτcgfOn the vertex by fsζτc=γfandgsζ=γggsζgf=γgγf=fsζgfτcgf
GMild ymax
Tight umax
Tight umax
τc<maxτchτcfζmin<ζh<ζhfζfhτcfhOn the vertex by fsζτc=γfandτc2=γh1hsζγf=fsζhfτchfτchf=12ζhfγh

Table 3.

Global optima of the constrained optimal servo control problem.

4. Design and evaluation of feasible constraints

The optimal solutions in Tables 2 and 3 are only true if the constraint set is such that a solution exists. Indeed, depending on the constraint set, the optimal solution may not have a feasible solution. Therefore, before applying the constrained optimal control formulation, either any given constraint set should first be screened quickly to determine its basic feasibility or a constraint set should be designed to be feasible.

4.1. Optimal regulatory control

Conditions for a feasible ymaxumax: For a given ymax, there is a minimum available umaxvalue below which the optimal control problem is not feasible. Figure 5 demonstrates the effects of the constraint specifications on the feasible region in ζτcspace. The constraint imposed by Eq. (5d) lies on a vertical line that shifts and bends rightward as umaxdecreases. The constraint given in Eq. (5c) shifts upward as umaxdecreases, whereas the constraint in Eq. (5b) shifts downward as ymaxdecreases. Note that the feasible region only exists when the constraint curve of umaxis below that of ymax.

Figure 5.

Effects of the constraint specifications y max , u max , and u max ′ on the feasible region.

As indicated in Figure 5, the feasible region reduces in size as umaxand ymaxdecrease, exhibiting continued reduction until ceasing to exist if the constraints are lower than some minimum allowable values. Therefore, there exists a tangent point ζtτctin ζτcspace, where the two constraint curves of umaxand ymaxmeet at a single point; this point equates to the smallest feasible umaxfor a given ymax(or the smallest feasible ymaxfor a given umax) for different specifications of umaxand ymax.

Let utmaxbe the smallest possible value of umax. utmaxcan be obtained when ζ=ζt. ζtcan be calculated by solving the following equation:

dhrζτc=dhrζγggr1=0E17

Once ζtis obtained, utmaxcan then be derived as follows:

umaxt=hrζtΔDE18

To sum up, if umaxutmax, the constraint set ymaxumaxis feasible; otherwise, it is infeasible.

Similarly, for a given umax, there is a minimum available ymaxvalue below which the optimal control problem is not feasible. Let ymaxtbe the smallest possible ymax. The values of ymaxtand ζtcan be obtained by solving the following system of equations simultaneously:

hrζγggr1=γhE19
dhrζγggr1=0E20

To sum up, if ymaxymaxt, then the constraint set ymaxumaxis feasible; otherwise, it is infeasible.

Feasible umaxfor a given feasible set ymaxumax: For a feasible ymaxumax, it can be intuitively inferred from the geometrical analysis of the constraint curves that any positive umaxΔDwill result in a feasible region if there is no intersection between the constraint curves of τcgrζ=γgand hrζτc=γh. Furthermore, if the intersection, ζgh, exists, the constraint, umax, is feasible if ζghζgf. To evaluate the existence of an intersection between the two constraint curves τcgrζ=γgand hrζτc=γh, it is important to calculate τc, the value of τcwhere the two constraint curves are secant when ζ. τcof each constraint curve can be obtained by solving the following equations:

τcgrζ=γgζE21a
hrζτc=γhζE21b

Because limζgrζ=1and limζhrζτc=1/τc=γh, τcof the two constraint curves are as follows:

τcg=γgE22a
τch=1γhE22b

A vertex ζghexists when τcgτch. Therefore,

γg1γhE23

which yields

ymaxKpΔDumaxΔD1E24

Overall, for a given feasible ymaxumax, umaxis feasible under either of the following conditions: (1) Eq. (24) is not satisfied or (2) Eq. (24) is satisfied and ζghζgf. Otherwise, umaxis not feasible and should be increased until one of the conditions is satisfied. Note that if ymaxKpΔDumaxΔD<1, no vertex point is formed by γg=τcgrζand γh=hrζτc, i.e., case D does not exist. In such a situation, for the purpose of evaluating the conditions presented in Table 2, any extremely large value can be assigned to ζgh.

Figure 6 illustrates a procedure applied to design a feasible constraint set ymaxumaxumaxand test its feasibility.

Figure 6.

Procedure to design and test a feasible constraint set for optimal regulatory PI control of integrating system.

4.2. Optimal servo control

It is clear from their approaching values of yt,ut,utas tthat for the constraints by ymaxand umaxto be feasible, they must be greater than ΔYspand ΔYsp/K, respectively, whereas the constraint by umaxcan be set to any nonnegative value. Figure 7 illustrates how the three constraint specifications affect the feasible region. The constraint γggsζvertically splits the region into two, while the constraints imposed by τc2γhhsζand γffsζτchave a similar shape in ζτcspace. It shows that for any feasible constraint set ymaxumaxumax, the feasible region is bounded below but unbounded in the upper side. This means that a decrease in ymax,umaxand umaxwill narrow down the feasible region delimited by the three constraints, but the feasible region will always exist. Moreover, the shape of the three constraints indicates the feasible region is always convex.

Figure 7.

Effects of the constraint specifications y max , u max , and u max ′ on the feasible region.

5. Closed-loop performance

5.1. Optimal regulatory control

Consider the following integrating process as

Gps=1sE25

Table 4 presents the examples of the seven possible aforementioned cases, as based on various constraint specifications. Simulations are carried out for weighting factors ωy=ωu=0.5.

ExampleCaseConstraint specificationPI parameter
ymaxumaxumaxKCτI
1A0.702.702.701.411.41
2B0.702.701.111.101.10
3C0.362.702.701.931.51
4D0.2852.702.102.100.69
5E0.701.1052.701.962.95
6F0.301.202.702.180.98
7G0.701.201.371.371.56

Table 4.

Constraint requirements and corresponding optimal PI parameters for regulatory control example.

Figure 8 presents the resulting process variable, yt, controller output, ut, and its rate of change, ut, for the seven examples. As can be seen from the figure, the PI controller designed by the proposed method not only yields optimal control performance, but also strictly satisfies the respective ymax, umax, and umaxconstraint requirements.

Figure 8.

Time responses of the system for cases A to G: regulatory system.

5.2. Optimal servo control

Consider the following integrating process

Gps=10sE26

Table 5 lists the examples of the seven possible cases based on various constraint specifications. Simulations are carried out for weighting factors, ωy=ωu=0.5. Figure 9 presents the time responses by the proposed optimal PI controller for the seven examples. As seen in the responses, the resulting optimal PI controllers not only provide the stable and optimal closed-loop responses, but also satisfy the ymax, umax, and umaxrequirements, strictly.

ExampleCaseConstraint specificationPI parameter
ymaxumaxumaxKCτI
1A1.20.51.51.411.414
2B1.20.50.50.791.58
3C1.030.51.51.521.46
4D1.030.51.031.511.47
5E1.20.31.51.532.42
6F1.030.411.51.381.61
7G1.20.41.172.151.83

Table 5.

Constraint requirements and corresponding optimal PI parameters for servo control example.

Figure 9.

Time responses of the system for cases A to G: servo system.

6. Conclusions

A novel analytical design approach is introduced for optimal regulatory and servo PI control subjected to operational constraints and examined to integrating processes. Owing to incisive parameterization, a complex constrained optimal control problem can be reformulated and converted to a simple algebraic form in the new design parameter ζτcspace, which allows finding the conditions and locations for the global optima by graphical analysis without having to rely on the numerical or black-box optimization effort. The proposed closed-form solution of the constrained optimal controller establishes a direct relationship between the control and plant parameters by which the optimal PI parameters can be obtained in an easy and quick manner. This approach also provides the following useful insights into how the control parameters affect the plant and how a feasible constraint set can be designed and checked in the constrained optimal control.

Acknowledgments

This research was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2015R1D1A3A01015621) and by the Priority Research Centers Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (2014R1A6A1031189).

Conflict of interest

The authors confirm there are no conflicts of interest.

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Rodrigue Tchamna and Moonyong Lee (September 12th 2018). Constraint Handling Optimal PI Controller Design for Integrating Processes: Optimization-Based Approach for Analytical Design, PID Control for Industrial Processes, Mohammad Shamsuzzoha, IntechOpen, DOI: 10.5772/intechopen.74301. Available from:

chapter statistics

201total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

Related Content

This Book

Next chapter

Decoupling Control and Soft Sensor Design for an Experimental Platform

By Thamiles Rodrigues de Melo, Nathália Arthur Brunet Monteiro, Danilo Pequeno, Jaidilson Jó da Silva and José Sérgio da Rocha Neto

Related Book

First chapter

Wavelet PID and Wavenet PID: Theory and Applications

By José Alberto Cruz Tolentino, Alejandro Jarillo Silva, Luis Enrique Ramos Velasco and Omar Arturo Domínguez Ramírez

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More about us