Globally convergent homotopy method for designing piecewise linear deterministic contractual function

In this paper, to design a piecewise linear contractual function, we 
consider to solve the single-level nonconvex programming with 
integral operator which is equivalent to the principal-agent bilevel 
programming model with continuous distribution. A modified 
constraint shifting homotopy method for solving the 
Karush-Kuhn-Tucker system of the discrete nonconvex programming is 
proposed and the global convergence from any initial point in 
shifted feasible set is proven under some mild conditions. A simple 
homotopy path tracing algorithm is given and is implemented in 
Matlab. For some typical risk averse utility functions and the 
typical distribution functions which simultaneously satisfy monotone 
likelihood ratio condition and convexity of the distribution 
function condition, some numerical tests to design the piecewise 
linear contract are done by our homotopy method as well as by using 
 fmincon in Matlab, LOQO and MINOS and, as a comparison, the 
piecewise constant contracts are also designed by solving the 
single-level nonconvex programming which is equivalent to the 
principal-agent bilevel programming model with corresponding 
discrete distributions. Numerical tests show that: to design a 
piecewise linear contract, which is much better than a piecewise 
constant contract, it needs only to solve a much lower dimensional 
optimization problem and hence needs much less computing time. 
Numerical experiences also show that the modified constraint 
shifting homotopy method is feasible and robust.

1. Introduction. In this paper, the principal-agent model with hidden action which was proposed by Mirrless [1] is considered. This is a bilevel programming problem, finite dimensional in case of discrete distribution or infinite dimensional in case of continuous distribution, in which the upper level object is that principal maximizes its expected utility through designing an incentive feasible compensation contract, and the constraints are the agent's individual rational constraint and incentive compatibility constraint, which is the lower level optimization problem, i.e., the agent privately chooses a level of effort to maximize its own expected utility.
To transform it into a single level programming problem, an effective and conventional method is first-order approach (FOA) which substitutes the lower level problem by its Karush-Kuhn-Tucker (KKT) condition. Mirrlees [1,2] and Rogerson [3] have proved the validity of FOA in the context of the assumptions that the distribution functions simultaneously satisfy monotone likelihood ratio condition(MLRC) and convexity of the distribution function condition(CDFC). Some distribution functions which simultaneously satisfy MLRC and CDFC have been constructed in [3,4,5]. For more references about first-order approach for solving principal-agent problems see, e.g., [6,7,8,9,10,11].
After transforming the bilevel programming into a single level programming in the context of satisfying MLRC and CDFC, one can try to solve it by the existing algorithms for nonlinear programming, with suitable design of the contract in the case of continuous distribution. However, generally speaking, the single-level programming transformed from bilevel programming by FOA is a nonconvex programming and it is difficult to find its global optimal solution. Even if one wants to find a local optimal solution or a KKT point of it, it needs some strong assumptions to get global convergence. In order to avoid the difficulty of solving nonconvex programming, Myerson [12] proposed the generalized principal-agent model in the case of discrete distribution, which is a linear programming problem by taking the probability distribution of each possible compensation, output and private information triplet as choice variable. Prescott [13,14] proposed the Dantzig-Wolfe decomposition algorithm for solving the generalized principal-agent model. A possible problem one may meet when he uses the excellent linear programming approach is: if the grid is sufficient fine, the size of the linear programming will be too big to solve. To avoid dealing with linear programming with huge size, Su and Judd [15] proposed a mathematical programming with equilibrium constraints (MPEC) approach and a hybrid procedure that combined MPEC and the linear programming approach for solving the generalized principal-agent model, in which they again met the difficulty of nonconvex programming and tried to find a solution as good as possible by using multi-initial points. In addition, solving generalized principal-agent model can only obtain the optimal probability distribution of each possible compensation, output and action triplet and not a specific contract.
To design a specific deterministic contract, we will try to give a globally convergent solution method for solving the equivalent nonconvex infinite-dimensional single-level nonlinear programming transformed from the bilevel programming of principal-agent model in continuous distribution case.
As an important globally convergent method, the homotopy method has been paid much attention for solving nonlinear problems since the 1970's, see, e.g., [16,17,18]. For solving nonconvex nonlinear programming with inequality constraint, in [19] and [20], a combined homotopy interior point method was proposed, and the existence and global convergence of a smooth homotopy path from a given interior point to a KKT point was proved under the normal cone condition as well as some commonly used constraint qualifications. To weaken the normal cone condition and permit non-interior initial point, Yu and Shang [21] constructed a constraint shifting homotopy. Lin, Li and Yu [22] generalized the combined homotopy method to nonconvex programming with both equality and inequality constraints. Yang, Yu and Xu [23] constructed a modified homotopy which is more convenient to use and the global convergence was proved under weaker conditions than the one given in [22].
In this paper, we consider to design a piecewise linear contractual function instead of designing general continuous contractual functions or designing discrete contractual functions, which in fact are piecewise constant contractual functions. To do so, we will solve the single-level nonconvex programming with integral operator which is equivalent to the principal-agent model in the space of piecewise linear functions. A globally convergent modified constraint shifting homotopy method for computing the solution to the KKT systems of the equivalent single level nonconvex programming with piecewise linear contractual function are proposed.
The remainder of this paper is organized as follows. In section 2, the infinite dimensional single-level nonconvex programming which is equivalent to the principalagent model with continuous distribution and its finite dimensional version via designing a piecewise linear contractual function are introduced. Nonconvexity of them are also discussed. In section 3, a modified constraint shifting homotopy method for solving KKT systems of the single level nonconvex programming with piecewise linear contractual function is proposed and the existence and convergence of a smooth path from any initial point in shifted feasible set to the solution to KKT systems is proved. In section 4, a simple homotopy path tracing algorithm is listed and some numerical tests to design the piecewise linear contract are done by our homotopy method as well as by using fmincon in Matlab, LOQO and MINOS and, as a comparison, the piecewise constant contracts are also designed by solving the single-level nonconvex programming which is equivalent to the principal-agent bilevel programming model with corresponding discrete distributions. Numerical tests show that: to design a piecewise linear contract, which is much better than a piecewise constant contract, it needs only to solve a lower dimensional optimization problem and hence needs much less computing time. Numerical experiences also show that the modified constraint shifting homotopy method is robust and effective. The conclusions are made in section 5.

2.1.
The principal-agent model. In the paper, we consider the following principal-agent model with hidden action in continuous distribution case which is a bilevel programming problem: where are the expected utilities of the principal and the agent, u(·) is the utility function of the principal, v(·) is the utility function of the agent, s(x) is the contractual function, i.e., the wage paid to the agent for the output x, (1b) is the individual rational constraint in which V 0 is the minimum expected utility of the agent, (1c) is the incentive compatibility constraint. The principal observes only a realized output x ∈ X = [x, x], the set of possible output, and the agent chooses an action a ∈ A = [a, a] on the basis of the agreed payment schedule s(x). f (x, a) is the continuous density function of output x for the action a undertaken by the agent, which is assumed sufficiently continuous differentiable in x and a. Let F (x, a) denote the distribution function of f (x, a).
If the contract s(x) is a piecewise constant on x and f (x, a) is a piecewise constant on x when the action a is taken by the agent, the model (1) is equivalent to the principal-agent model in discrete case, see [3].
As in the literatures, suppose that the principal is risk neutral with utility u(x − s(x)) = x − s(x), and the agent is risk averse with a separable Von Neumann- ) < 0 and c (a) ≥ 0. And suppose that the following two properties hold: (I) f (x, a) satisfies MLRC, i.e., fa(x,a) f (x,a) is nondecreasing in x for every a. MLRC is used to insure that s(x) is nondecreasing in outcome x. (II) f (x, a) satisfies CDFC, i.e., F aa (x, a) ≥ 0. CDFC is used to insure that the expected utility of the agent is concave in a. Then, by FOA, we have the following equivalent relaxed pareto-optimization programming of the principal-agent model with continuous distribution: in which constraint (1c) is replaced with (4c). When the contract s(x) is a piecewise constant on x and f (x, a) is a piecewise constant on x when the action a is taken by the agent, the MLRC and CDFC are equivalent to the corresponding properties of the discrete distributions and the single-level programming (4) is also equivalent to the corresponding relaxed paretooptimization programming in discrete case, see [3].
Generally speaking, the equivalent relaxed pareto-optimization programming (4) is a nonconvex nonlinear programming. Firstly, the equality constraint (4c) is generally nonlinear, consequently problem (4) is nonconvex. For example, let utility function of the risk averse agent be v(s(x)) = − exp(−s(x)) and the distribution function be F (x, a) = x * exp(a(x − 1)) which is constructed in [4], then Obviously, the equality constraint function V a (s(x), a) is nonlinear.
Secondly, taking second order partial derivative with s(x) of the inequality constraint function V (s(x), a) in (4), we have Since v < 0 and f (x, a) is a density function, we get V ss (s(x), a) ≤ 0. Therefore, V (s(x), a) is concave in s(x). Although the inequality constraint function V (s(x), a) is concave in both s(x) and a implied by CDFC, V (s(x), a) is generally a nonconcave function. For example, let G(x, y) = −xy + 1, although G(x, y) is not only concave but also convex in x and y, it is neither concave nor convex. Since the Hessian matrix of G(x, y) is By simple computation, the eigenvalues of the Hessian matrix ∇ 2 G are −1 and 1. Then, the Hessian matrix ∇ 2 G is indefinite, i.e., G(x, y) is neither concave nor convex.
Then, the expected utilities (2) and (3) of the principal-agent model with hidden action can be reformulated as follows: Therefore, the relaxed pareto-optimization programming (4) of the principalagent model with piecewise linear deterministic contractual function can be written as follows:
As it's known that it is difficult to obtain the global solution of nonconvex nonlinear programming, we can only obtain the solution to its KKT systems. Therefore, in order to solve problem (9), we need to solve the following KKT systems: where z ∈ R m , ∇ = ( ∂ ∂p , ∂ ∂q , ∂ ∂a ) T . Throughout the paper, let R n , R n + , R n ++ denote n-dimension Euclidean space, nonnegative orthant and positive orthant of R n , respectively. The following parametrized Sard's lemma (see Theorem 2.1., [18]) will be used in this paper. Let B ⊂ R n be an open set and let ϕ : B → R p be a C α (α > max{0, n − p}) mapping; we say that π ∈ R p is a regular value of ϕ if Lemma 2.1. Let B ⊂ R n , D ⊂ R m be open sets, and let ϕ : B × D → R k be a C α mapping, where α > max{0, m − k}. If 0 ∈ R k is a regular value of ϕ, then for almost all b ∈ B, 0 is a regular value of ϕ b = ϕ(b, ·).
(A5) (The normal cone condition) ∀θ ∈ ∂Ω(1), the normal cone of Ω(1) only meets ∂Ω(1) at θ, i.e., ∀θ ∈ ∂Ω(1), Remark 1. In comparison with [22], the normal cone condition (A5) only requires that the constraint shifting set Ω(t) holds as t = 1 but not the whole constraint set Ω. Therefore, it is more weaker than that of [22]. In comparison with [23], the construction of the shifted constraint functions in homotopy (11) is more optional and our method can solve the more general problems, as is shown in the following example that some problems can be solved by our method while which can't be solved by the method in [23].
are two discs. The nonconvex feasible set Ω which is drawn in the Figure 1(a), 2(a), 3(a) and 4(a) (the area filled with gray) doesn't satisfy the normal cone condition, since the normal line drawn in the figure intersects the set Ω itself.
Since the normal line intersects the set Ω(1) itself as is shown in the Figure 1(d), the shifting feasible set Ω(1) doesn't satisfy normal cone condition. Moreover, as is shown in the Figure 1(c), for t = 0.31 the shifting feasible set Ω(t) doesn't satisfy the positive linearly independent condition at the point, at which the smaller disc D 2 (t) is internally tangent to the bigger disc D 1 (t) and the two normal lines are in the opposite direction. Case 2. When ρ = 1 and τ = 2e 2 = 2 * (0, 1) T , the shifted constraint function As parameter t changes from t = 1 to t = 0, the shifting feasible set (13) will deform from Ω(1) to Ω(0) = Ω, in which the bigger disc D 1 (t) = D 1 will keep fixed, while the smaller disc D 2 (t) = {x : −ḡ 2 (x, t) < 0} will gradually expand from an empty set D 2 (1) to D 2 . For more details about the shifting process see Figure 2 (a)-(d) for t = 0, 0.18, 0.4 and 1, respectively. As is shown in the Figure 2(d), for t = 1 the shifting feasible set Ω(t) is the whole disc D 1 (1) itself which is a convex set, naturally, it satisfies the normal cone condition. However, as is shown in the Figure 2(b), for t = 0.18 the shifting feasible set Ω(t) doesn't satisfy the positive linearly independent condition at the point, at which the smaller disc D 2 (t) is internally tangent to the bigger disc D 1 (t) and the two normal lines are in the opposite direction. In comparison with the shifted constraint functionḡ(x, t) given in [23], we can construct a shifted constraint function of the formg(x, t) given in the modified homotopy (11) to satisfy both conditions (A4) and (A5). Two typical cases are given below to show its feasibility. Case A. The shifted constraint functiong(x, t) is constructed as In the shifting process, the shifting feasible set (14) will deform from Ω (1) to Ω(0) = Ω as parameter t changes from t = 1 to t = 0, in which the bigger disc D 1 (t) = {x :g 1 (x, t) ≤ 0} will gradually expand up to D 1 while the smaller disc D 2 (t) = D 2 will keep fixed. For more details about the shifting process see Figure 3 (a)-(d) for t=0, 0.2, 0.6 and 1, respectively.  Figure 3. The process of constraint shifting generated by (14) for t=0,0.2,0.6,1.
As is shown in the Figure 3(d), for t = 1 any normal line doesn't intersect the shifting feasible set Ω(t) itself, then the shifting feasible set Ω(1) satisfies normal cone condition. Moreover, in the shifting process as is shown in the Figure 3, the shifting feasible set Ω(t) always satisfies the positive linearly independent condition when t changes from t = 1 to t = 0. Case B. The shifted constraint functiong(x, t) is constructed as When parameter t changes from t = 1 to t = 0, the shifting feasible set Ω(t) = D 1 (t) \ D 2 (t) = {x :g(x, t) ≤ 0} generated by (15) will deform from Ω(1) to Ω(0) = Ω, in which the disc D 1 (t) = {x : (x 1 + t) 2 + x 2 2 − 9 ≤ 0} and D 2 (t) = {x : −g 2 (x, t) < 0} will keep the shape unchanged and move towards each other. For more details about the shifting process see Figure 4    As is shown in the Figure 4(d), for t = 1 the shifting feasible set Ω(t) is the whole disc D 1 (1) itself which is a convex set, naturally, it satisfies normal cone condition. Moreover, in the shifting process as is shown in the Figure 4 (a)-(d), Ω(t) always satisfies the positive linearly independent condition as t changes from t = 1 to t = 0.
For the modified homotopy (11), when t = 0, the homotopy equation H(w, 0) = 0 turns to the KKT systems (10) of problem (9); when t = 1, the homotopy equation H(w, 1) = 0 has a unique simple solution, which can be proved by the following lemma.
Therefore, we obtain the result.
If the case (ii) happens, there exists a sequence of points {(θ k , y k , z k , t k )} ⊂ Γ w 0 such that t k → t * ∈ [0, 1], θ k → θ * ∈ Ω(t * ) and (y k , z k ) → ∞ as k → ∞. From the first equation of (11), we have And only the following two subcases are possible: ( , then y * i = 0 for i / ∈ I 1 (θ * ) from the second equation of (11). Taking limits in (16) as k → ∞, we have which contradicts with condition (A5). If ((1 − t k )y k , z k ) → ∞, the discussion is the same as case (II).
As a conclusion, case (v) is the only possible case. That is Γ w 0 must terminate in or approach to the hyperplane at t * = 0 and w * = (θ * , y * , z * ) is a solution to KKT systems (10) of problem (9).
Remark 2. In Theorem 3.2, the global convergence only means that the proposed method, i.e., the modified constraint shifting homotopy method, converges globally for solving the KKT systems (10).

Implementation of the algorithm and numerical examples.
For the modified constraint shifting homotopy (11), the numerical path-tracing of the homotopy path Γ w 0 can be implemented by the following procedure, in which predictorcorrector procedure contains three steps: the first predictor step, the midway predictor step and the corrector step. The first predictor step is taken by computing the tangent direction, and the midway predictor steps are taken by using secant directions. The corrector steps are taken by Newton iterations for solving an augmented system. Some detailed discussions on the predictor-corrector algorithms and the convergence can be seen, e.g., [25,26] and etc. Detailed description of the implementation on homotopy method is formulated as follows.
Step 2. The first predictor step.
Step 5. The midway predictor step.
Step 6. The end game.

ZHICHUAN ZHU, BO YU AND LI YANG
Compute the Newton step d end by solving the equation If H(w k+1,j , 0) ∞ ≤ ε 2 and w k+1,j is feasible, then terminate the algorithm. w = w k+1,j is the computed solution to the KKT systems (10). Therefore, a = a k+1,j is the optimal action and s k+1,j i (x) = p k+1,j i x+q k+1,j i , i = 1, 2, . . . , m is the piecewise linear contractual function.
Algorithm 4.1 cannot only be applied to design piecewise linear contractual function but also to design piecewise constant contractual function. We have implemented it by Matlab language and done some numerical tests for some typical utility functions and distribution functions satisfying MLRC and CDFC. The computations are performed on a computer running the software Matlab R2007b on Microsoft Windows XP Professional with Intel(R) 1.83GHz processor and 1024 megabytes of memory.
To show the advantage of designing piecewise linear contractual function over designing piecewise constant contractual function, we have made comparison between numerical results for both models. To show the feasibility, robustness and effectiveness of Algorithm 4.1, we have also computed the test problems by Matlab minimization routine fmincon, LOQO 6.01 Student Version and MINOS 5.5 Student Version.
For a sake of simpleness, in the following numerical examples, the shifted constraint functions are constructed asg i (θ, t) = g i (θ) − t 2 τ andh j (θ, t) = h j (θ) − t 2 h j (θ 0 ), where the constant τ is chosen as τ = 0 when g i (θ 0 ) < 0 and τ = max{g i (θ 0 )} + 10 when g i (θ 0 ) ≮ 0. Of course, we can also construct special shifted constraint functions based on the feature of the usual risk averse utility functions. For example, for the typical risk averse utility function v(s(x)) = − exp(−γ * s(x)), where γ > 0 is the coefficient of risk aversion, if the solution to the principalagent model with γ = γ 1 is easy to compute (or has been computed), then for the principal-agent model with γ = γ 2 , we can construct the shifted constraint function asg(θ, t) = x x − exp(−γ(t) * s(x))f (x, a)dx − a, where γ(t) = tγ 1 + (1 − t)γ 2 and t is homotopy parameter.
In modified homotopy (11) and Algorithm 4.1, the parameters are set as follows: η i = 10, ξ = randn(2m + 1, 1), a 2m + 1 vector with random entries drawn from a normal distribution with mean zero and standard deviation one,N =5, λ 0 =0.7, 2, ε 1 = 10 −2 and ε 2 = 10 −6 . In the following numerical tests, designing piecewise linear contractual function is made comparison with designing piecewise constant contractual function as the following three aspects: (1) when the discrete segment is same, the optimal objective value of designing piecewise linear contractual function is made comparison with designing piecewise constant contractual function; (2) when the computing time is almost same, how many segments need to be discretized is discussed in case of designing piecewise constant contractual function; (3) when the optimal objective value is almost same, how many segments need to be discretized is also discussed in case of designing piecewise constant contractual function.
In the following tables, the computing methods are listed in column 1, the contract is listed in column 2, the optimal objective value U * is listed in column 3, the computing time denoted by CPU is listed in column 4, the iteration steps denoted by IT is listed in column 5. Let MCSH denote the modified constraint shifting homotopy, PL-m and PC-m denote piecewise linear contractual function and piecewise constant contract with the segment number m, respectively, e.g., PL-2 and PC-2 denote piecewise linear contractual function and piecewise constant contract with the segment number m=2.
Since LOQO 6.01 student version and MINOS 5.5 Student Version are limited to 300 variables and 300 constraints and objectives, and iteration steps are limited to 500, when the variables are more than 300 or iteration steps exceed 500, LOQO 6.01 student version and MINOS 5.5 Student Version are invalid and it will be indicated by " "in the following tables. For the MINOS 5.5 Student Version, if some functions such as pow(·) and exp(·) can't be evaluated and lead to the solution aborted, it will be indicated by " "in the following tables.
Since the limitations of fmincon, if the number of iterations exceeds maximum number 99 of iterations allowed and the measure of first-order optimality condition is too large or infinity, fmincon will be terminated and it will be indicated by "--". If a method terminates in a finite iteration, but it doesn't satisfy the final precision required and the result isn't a solution to KKT systems, it will be indicated by "!". If a method goes into an endless loop or fails, it will be indicated by dash "-". If a method can't show the computing time, it will be indicated by " ".
For the method MCSH and the initial points θ 0 1 , θ 0 2 and θ 0 3 , when the contract is piecewise linear function, the numerical solution is θ * =(1.688,0.815,0.073,0.013, 0.384). Therefore, the action is a * =1.688 and the contract is s(x) = 0.815x + 0.013, x 1 = 0 < x ≤ 0.50 0.073x + 0.384, 0.50 ≤ x < x 3 = 1 which is drawn in the Figure 5.    From Table 4-6, when the contract is piecewise linear function, in comparison with fmincon, since the optimal objective values of fmincon are 0.3355 and 0.3356 which are different and smaller than 0.3357 for the initial points θ 0 2 and θ 0 3 , MCSH, LOQO 6.01 and MINOS 5.5 are more robust and effective. When the contract is piecewise constant and the segment number is 2 which is equal to that of designing piecewise linear contractual function, the optimal objective values is only 0.0652 for the four methods and MINOS 5.5 fails for the second initial pointsθ 0 2 ; When the computing time is almost same, designing piecewise constant contract  needs to be discretized as 70 segments and the optimal objective values is 0.3285 for the four methods except MINOS 5.5; When the optimal objective value of designing piecewise constant contract is almost equal to that of designing piecewise linear contractual function, the piecewise constant contract needs to be discretized as 6000 segments, the computing time is more than 600 seconds and the results are concluded as follows: the optimal objective value of MCSH is 0.3356, the optimal objective value of fmincon is 0.3354, 0.3331 and 0.3351 for the initial pointsθ 0 1 ,θ 0 2 andθ 0 3 respectively, LOQO 6.01 and MINOS 5.5 are invalid since the variables exceed 300. Therefore, designing piecewise linear contractual function is much better than designing piecewise constant contract.
For the method MCSH and the initial points θ 0 1 , θ 0 2 and θ 0 3 , when the contract is piecewise linear function, the numerical solution is θ * =(0.906,0.368,0.494,0.0089,-0.054). Therefore, the action is a * =0.906 and the contract is which is drawn in the Figure 6. Example 4. For distribution function F (x, a) = x · exp(a(x − 1)) which is constructed in [4], x ∈ (0, 1), a ∈ A = R + , c(a) = 0.06a and the utility function of agent v(s(x)) = − exp(−s(x)), suppose the piecewise linear contract function is  Table 7-9. From Table 7-9, when the contract is piecewise linear function, the measure of first-order optimality condition of fmincon is infinity for θ 0 1 and θ 0 2 , and the optimal objective value of fmincon is 0.0239 for θ 0 3 ; MINOS 5.5 is solution aborted for θ 0 1 and θ 0 2 ; However, the optimal objective values of MCSH and LOQO 6.01 are  0.0240 for the three initial points, which show that MCSH and LOQO 6.01 are much better. When the contract is piecewise constant, for the segment number 4 which is equal to that of designing piecewise linear contractual function, the optimal objective values is only -0.1078 for the four methods which is much smaller than 0.0240; When the computing time is almost same, the piecewise constant contract needs to be discretized as 85, 240 and 70 segments for the three initial points and the optimal objective values are still much smaller than 0.0240; When the optimal objective value of designing piecewise constant contract is almost equal to 0.0240, for the three initial points, the piecewise constant contract needs to be discretized as 6000 segments, the computing time needs at least 350 seconds and the results are concluded as follows: the optimal objective value of MCSH is 0.0239, fmincon can't obtain a solution to KKT systems when the procedure terminates, LOQO 6.01 and MINOS 5.5 are invalid since the variables exceed 300. Therefore, the numerical results show that designing piecewise linear contractual function is better than designing piecewise constant contract, and MCSH is more effective and feasible. For the method MCSH and the initial points θ 0 1 , θ 0 2 and θ 0 3 , when the contract is piecewise linear function, the numerical solution is θ * =(0.62,1.27,0.86,0.66,0.53, 0.004,0.10,0.21,0.30). Therefore, the action is a * = 0.62 and the contract is  180 seconds and the results are concluded as follows: the optimal objective value of method MCSH is 0.2447, fmincon can't obtain a solution to KKT systems when the procedure terminates, LOQO 6.01 and MINOS 5.5 are invalid since the variables exceed 300. Therefore, the numerical results shows that the piecewise linear contract is better than piecewise constant contract, and the method MCSH is more robust and effective. For the method MCSH and the initial points θ 0 1 , θ 0 2 and θ 0 3 , when the contract is piecewise linear function, the numerical solution is θ * =(1.51,0.29,0.33,0.37,0.42,0.49, 0.57,0.68,0.82,0.02,0.02,0.006,-0.01,-0.05,-0.10,-0.18,-0.30). Therefore, the action a * =1.51 and the piecewise linear contractual function is 8 ≤ x < x 9 = 1 which is drawn in Figure 8.

5.
Conclusions. The main work of this paper is listed as follows: (1) to design a piecewise linear contractual function, we consider to solve the single-level infinitedimensional nonconvex programming with integral operator which is equivalent to the principal-agent bilevel programming model; (2) the globally convergent modified constraint shifting homotopy method is proposed in order to solve the KKT systems of the equivalent relaxed pareto-optimization programming with piecewise linear deterministic contractual function and the global convergence is obtained for any initial point chosen in shifting feasible set under some mild conditions; (3) in the case of typical risk averse utility function and typical distribution functions simultaneously satisfying monotone likelihood ratio condition and convex of distribution function condition, some numerical tests to design the piecewise linear contract are done by our homotopy method as well as by using fmincon in Matlab, LOQO and MINOS and, as a comparison, the piecewise constant contracts are also designed by solving the single-level nonconvex programming which is equivalent to the principal-agent bilevel programming model with corresponding discrete distributions.

ZHICHUAN ZHU, BO YU AND LI YANG
The numerical tests in Section 4 show that: 1) As is shown in the figures in Section 4, the piecewise linear contractual function is nondecreasing in output x. This coincides with the theoretical result (Prop. 1 (iii), [3]); 2) In comparison with designing piecewise constant contract, the advantage of designing a piecewise linear contract is that it needs only to solve a much lower dimensional optimization problem and needs much less computing time to get a similar optimal value; 3) For the problems that can be solved by Matlab minimization routine fmincon, LOQO and MINOS, our modified constraint shifting homotopy method can also solve them and get the same optimal objective values and solutions. The choice of the initial points of the modified constraint shifting homotopy method is flexible and convenient. So we can say that the modified constraint shifting homotopy method is feasible, effective and robust.