SSR: SemiSmooth Reformulation

Description

functions of the SemiSmooth Reformulation of the GNEP

Usage

funSSR(z, dimx, dimlam, grobj, arggrobj, constr, argconstr,  grconstr, arggrconstr, 
	compl, argcompl, dimmu, joint, argjoint, grjoint, arggrjoint, echo=FALSE)
jacSSR(z, dimx, dimlam, heobj, argheobj, constr, argconstr,  grconstr, arggrconstr, 
	heconstr, argheconstr, gcompla, gcomplb, argcompl, dimmu, joint, argjoint,
	grjoint, arggrjoint, hejoint, arghejoint, echo=FALSE)

Value

A vector for funSSR or a matrix for jacSSR.

Arguments

z: a numeric vector z containing \((x, lambda, mu)\) values.
dimx: a vector of dimension for \(x\).
dimlam: a vector of dimension for \(lambda\).
grobj: gradient of the objective function, see details.
arggrobj: a list of additional arguments of the objective gradient.
constr: constraint function, see details.
argconstr: a list of additional arguments of the constraint function.
grconstr: gradient of the constraint function, see details.
arggrconstr: a list of additional arguments of the constraint gradient.
compl: the complementarity function with (at least) two arguments: compl(a,b).
argcompl: list of possible additional arguments for compl.
dimmu: a vector of dimension for \(mu\).
joint: joint function, see details.
argjoint: a list of additional arguments of the joint function.
grjoint: gradient of the joint function, see details.
arggrjoint: a list of additional arguments of the joint gradient.
heobj: Hessian of the objective function, see details.
argheobj: a list of additional arguments of the objective Hessian.
heconstr: Hessian of the constraint function, see details.
argheconstr: a list of additional arguments of the constraint Hessian.
gcompla: derivative of the complementarity function w.r.t. the first argument.
gcomplb: derivative of the complementarity function w.r.t. the second argument.
hejoint: Hessian of the joint function, see details.
arghejoint: a list of additional arguments of the joint Hessian.
echo: a logical to show some traces.

Author

Christophe Dutang

Details

Compute the SemiSmooth Reformulation of the GNEP: the Generalized Nash equilibrium problem is defined by objective functions \(Obj\) with player variables \(x\) defined in dimx and may have player-dependent constraint functions \(g\) of dimension dimlam and/or a common shared joint function \(h\) of dimension dimmu, where the Lagrange multiplier are \(lambda\) and \(mu\), respectively, see F. Facchinei et al.(2009) where there is no joint function.

Arguments of the Phi function

The arguments which are functions must respect the following features

grobj: The gradient \(Grad Obj\) of an objective function \(Obj\) (to be minimized) must have 3 arguments for \(Grad Obj(z, playnum, ideriv)\): vector z, player number, derivative index , and optionnally additional arguments in arggrobj.
constr: The constraint function \(g\) must have 2 arguments: vector z, player number, such that \(g(z, playnum) <= 0\). Optionnally, \(g\) may have additional arguments in argconstr.
grconstr: The gradient of the constraint function \(g\) must have 3 arguments: vector z, player number, derivative index, and optionnally additional arguments in arggrconstr.
compl: It must have two arguments and optionnally additional arguments in argcompl. A typical example is the minimum function.
joint: The constraint function \(h\) must have 1 argument: vector z, such that \(h(z) <= 0\). Optionnally, \(h\) may have additional arguments in argjoint.
grjoint: The gradient of the constraint function \(h\) must have 2 arguments: vector z, derivative index, and optionnally additional arguments in arggrjoint.

Arguments of the Jacobian of Phi

The arguments which are functions must respect the following features

heobj: It must have 4 arguments: vector z, player number, two derivative indexes and optionnally additional arguments in argheobj.

heconstr

It must have 4 arguments: vector z, player number, two derivative indexes and optionnally additional arguments in argheconstr.

gcompla,gcomplb

It must have two arguments and optionnally additional arguments in argcompl.

hejoint

It must have 3 arguments: vector z, two derivative indexes and optionnally additional arguments in arghejoint.

See the example below.

References

F. Facchinei, A. Fischer and V. Piccialli (2009), Generalized Nash equilibrium problems and Newton methods, Math. Program.

Examples

Run this code


# (1) associated objective functions
#

dimx <- c(2, 2, 3)

#Gr_x_j O_i(x)
grfullob <- function(x, i, j)
{
	x <- x[1:7]	
	if(i == 1)
	{
		grad <- 3*(x - 1:7)^2
	}
	if(i == 2)
	{
		grad <- 1:7*(x - 1:7)^(0:6)
	}
	if(i == 3)
	{
		s <- x[5]^2 + x[6]^2 + x[7]^2 - 5	
		grad <- c(1, 0, 1, 0, 4*x[5]*s, 4*x[6]*s, 4*x[7]*s)
			
	}
	grad[j]	
}


#Gr_x_k Gr_x_j O_i(x)
hefullob <- function(x, i, j, k)
{
	x <- x[1:7]
	if(i == 1)
	{
		he <- diag( 6*(x - 1:7) )
	}
	if(i == 2)
	{
		he <- diag( c(0, 2, 6, 12, 20, 30, 42)*(x - 1:7)^c(0, 0:5) )
	}
	if(i == 3)
	{
		s <- x[5]^2 + x[6]^2 + x[7]^2	
		
		he <- rbind(rep(0, 7), rep(0, 7), rep(0, 7), rep(0, 7),
			c(0, 0, 0, 0, 4*s+8*x[5]^2, 8*x[5]*x[6], 8*x[5]*x[7]),
			c(0, 0, 0, 0, 8*x[5]*x[6], 4*s+8*x[6]^2, 8*x[6]*x[7]),
			c(0, 0, 0, 0,  8*x[5]*x[7], 8*x[6]*x[7], 4*s+8*x[7]^2) )
	}
	he[j,k]	
}



# (2) constraint linked functions
#

dimlam <- c(1, 2, 2)

#constraint function g_i(x)
g <- function(x, i)
{
	x <- x[1:7]
	if(i == 1)
		res <- sum( x^(1:7) ) -7
	if(i == 2)
		res <- c(sum(x) + prod(x) - 14, 20 - sum(x))
	if(i == 3)
		res <- c(sum(x^2) + 1, 100 - sum(x))
	res
}


#Gr_x_j g_i(x)
grfullg <- function(x, i, j)
{
	x <- x[1:7]	
	if(i == 1)
	{
		grad <- (1:7) * x ^ (0:6)
	}
	if(i == 2)
	{
		grad <- 1 + sapply(1:7, function(i) prod(x[-i]))
		grad <- cbind(grad, -1)
	}
	if(i == 3)
	{
		grad <- cbind(2*x, -1)
	}


	if(i == 1)
		res <- grad[j]	
	if(i != 1)
		res <- grad[j,]	
	as.numeric(res)
}



#Gr_x_k Gr_x_j g_i(x)
hefullg <- function(x, i, j, k)
{
	x <- x[1:7]
	if(i == 1)
	{
		he1 <- diag( c(0, 2, 6, 12, 20, 30, 42) * x ^ c(0, 0, 1:5) )
	}
	if(i == 2)
	{
		he1 <- matrix(0, 7, 7)
		he1[1, -1] <- sapply(2:7, function(i) prod(x[-c(1, i)]))
		he1[2, -2] <- sapply(c(1, 3:7), function(i) prod(x[-c(2, i)]))
		he1[3, -3] <- sapply(c(1:2, 4:7), function(i) prod(x[-c(3, i)]))
		he1[4, -4] <- sapply(c(1:3, 5:7), function(i) prod(x[-c(4, i)]))
		he1[5, -5] <- sapply(c(1:4, 6:7), function(i) prod(x[-c(5, i)]))
		he1[6, -6] <- sapply(c(1:5, 7:7), function(i) prod(x[-c(6, i)]))
		he1[7, -7] <- sapply(1:6, function(i) prod(x[-c(7, i)]))
						
						
		he2 <- matrix(0, 7, 7)
		
	}
	if(i == 3)
	{
		he1 <- diag(rep(2, 7))
		he2 <- matrix(0, 7, 7)
	}
	if(i != 1)
		return( c(he1[j, k], he2[j, k])	)
	else				
		return( he1[j, k] )
}


# (3) compute Phi
#

z <- rexp(sum(dimx) + sum(dimlam))

n <- sum(dimx)
m <- sum(dimlam)
x <- z[1:n]
lam <- z[(n+1):(n+m)]

resphi <- funSSR(z, dimx, dimlam, grobj=grfullob, constr=g, grconstr=grfullg, compl=phiFB)


check <- c(grfullob(x, 1, 1) + lam[1] * grfullg(x, 1, 1), 
	grfullob(x, 1, 2) + lam[1] * grfullg(x, 1, 2), 
	grfullob(x, 2, 3) + lam[2:3] %*% grfullg(x, 2, 3),
	grfullob(x, 2, 4) + lam[2:3] %*% grfullg(x, 2, 4), 	
	grfullob(x, 3, 5) + lam[4:5] %*% grfullg(x, 3, 5),
	grfullob(x, 3, 6) + lam[4:5] %*% grfullg(x, 3, 6),
	grfullob(x, 3, 7) + lam[4:5] %*% grfullg(x, 3, 7),
	phiFB( -g(x, 1), lam[1]), 
	phiFB( -g(x, 2)[1], lam[2]), 
	phiFB( -g(x, 2)[2], lam[3]), 
	phiFB( -g(x, 3)[1], lam[4]), 
	phiFB( -g(x, 3)[2], lam[5]))
	
	

#check
cat("\n\n________________________________________\n\n")

#part A
print(cbind(check, res=as.numeric(resphi))[1:n, ])
#part B
print(cbind(check, res=as.numeric(resphi))[(n+1):(n+m), ])
	
# (4) compute Jac Phi
#
	
resjacphi <- jacSSR(z, dimx, dimlam, heobj=hefullob, constr=g, grconstr=grfullg, 
	heconstr=hefullg, gcompla=GrAphiFB, gcomplb=GrBphiFB)

	
#check
cat("\n\n________________________________________\n\n")


cat("\n\npart A\n\n")	


checkA <- 
rbind(
c(hefullob(x, 1, 1, 1) + lam[1]*hefullg(x, 1, 1, 1), 
hefullob(x, 1, 1, 2) + lam[1]*hefullg(x, 1, 1, 2),
hefullob(x, 1, 1, 3) + lam[1]*hefullg(x, 1, 1, 3),
hefullob(x, 1, 1, 4) + lam[1]*hefullg(x, 1, 1, 4),
hefullob(x, 1, 1, 5) + lam[1]*hefullg(x, 1, 1, 5),
hefullob(x, 1, 1, 6) + lam[1]*hefullg(x, 1, 1, 6),
hefullob(x, 1, 1, 7) + lam[1]*hefullg(x, 1, 1, 7)
),
c(hefullob(x, 1, 2, 1) + lam[1]*hefullg(x, 1, 2, 1), 
hefullob(x, 1, 2, 2) + lam[1]*hefullg(x, 1, 2, 2),
hefullob(x, 1, 2, 3) + lam[1]*hefullg(x, 1, 2, 3),
hefullob(x, 1, 2, 4) + lam[1]*hefullg(x, 1, 2, 4),
hefullob(x, 1, 2, 5) + lam[1]*hefullg(x, 1, 2, 5),
hefullob(x, 1, 2, 6) + lam[1]*hefullg(x, 1, 2, 6),
hefullob(x, 1, 2, 7) + lam[1]*hefullg(x, 1, 2, 7)
),
c(hefullob(x, 2, 3, 1) + lam[2:3] %*% hefullg(x, 2, 3, 1), 
hefullob(x, 2, 3, 2) + lam[2:3] %*% hefullg(x, 2, 3, 2),
hefullob(x, 2, 3, 3) + lam[2:3] %*% hefullg(x, 2, 3, 3),
hefullob(x, 2, 3, 4) + lam[2:3] %*% hefullg(x, 2, 3, 4),
hefullob(x, 2, 3, 5) + lam[2:3] %*% hefullg(x, 2, 3, 5),
hefullob(x, 2, 3, 6) + lam[2:3] %*% hefullg(x, 2, 3, 6),
hefullob(x, 2, 3, 7) + lam[2:3] %*% hefullg(x, 2, 3, 7)
),
c(hefullob(x, 2, 4, 1) + lam[2:3] %*% hefullg(x, 2, 4, 1), 
hefullob(x, 2, 4, 2) + lam[2:3] %*% hefullg(x, 2, 4, 2), 
hefullob(x, 2, 4, 3) + lam[2:3] %*% hefullg(x, 2, 4, 3), 
hefullob(x, 2, 4, 4) + lam[2:3] %*% hefullg(x, 2, 4, 4), 
hefullob(x, 2, 4, 5) + lam[2:3] %*% hefullg(x, 2, 4, 5), 
hefullob(x, 2, 4, 6) + lam[2:3] %*% hefullg(x, 2, 4, 6), 
hefullob(x, 2, 4, 7) + lam[2:3] %*% hefullg(x, 2, 4, 7)
),
c(hefullob(x, 3, 5, 1) + lam[4:5] %*% hefullg(x, 3, 5, 1),  
hefullob(x, 3, 5, 2) + lam[4:5] %*% hefullg(x, 3, 5, 2),  
hefullob(x, 3, 5, 3) + lam[4:5] %*% hefullg(x, 3, 5, 3),  
hefullob(x, 3, 5, 4) + lam[4:5] %*% hefullg(x, 3, 5, 4),  
hefullob(x, 3, 5, 5) + lam[4:5] %*% hefullg(x, 3, 5, 5),  
hefullob(x, 3, 5, 6) + lam[4:5] %*% hefullg(x, 3, 5, 6),  
hefullob(x, 3, 5, 7) + lam[4:5] %*% hefullg(x, 3, 5, 7)
),
c(hefullob(x, 3, 6, 1) + lam[4:5] %*% hefullg(x, 3, 6, 1),   
hefullob(x, 3, 6, 2) + lam[4:5] %*% hefullg(x, 3, 6, 2),  
hefullob(x, 3, 6, 3) + lam[4:5] %*% hefullg(x, 3, 6, 3),  
hefullob(x, 3, 6, 4) + lam[4:5] %*% hefullg(x, 3, 6, 4),  
hefullob(x, 3, 6, 5) + lam[4:5] %*% hefullg(x, 3, 6, 5),  
hefullob(x, 3, 6, 6) + lam[4:5] %*% hefullg(x, 3, 6, 6),  
hefullob(x, 3, 6, 7) + lam[4:5] %*% hefullg(x, 3, 6, 7)
),
c(hefullob(x, 3, 7, 1) + lam[4:5] %*% hefullg(x, 3, 7, 1),   
hefullob(x, 3, 7, 2) + lam[4:5] %*% hefullg(x, 3, 7, 2),  
hefullob(x, 3, 7, 3) + lam[4:5] %*% hefullg(x, 3, 7, 3),  
hefullob(x, 3, 7, 4) + lam[4:5] %*% hefullg(x, 3, 7, 4),  
hefullob(x, 3, 7, 5) + lam[4:5] %*% hefullg(x, 3, 7, 5),  
hefullob(x, 3, 7, 6) + lam[4:5] %*% hefullg(x, 3, 7, 6),  
hefullob(x, 3, 7, 7) + lam[4:5] %*% hefullg(x, 3, 7, 7)
)
)


print(resjacphi[1:n, 1:n] - checkA)


cat("\n\n________________________________________\n\n")


cat("\n\npart B\n\n")	


checkB <- 
rbind(
cbind(c(grfullg(x, 1, 1), grfullg(x, 1, 2)), c(0, 0), c(0, 0), c(0, 0), c(0, 0)),
cbind(c(0, 0), rbind(grfullg(x, 2, 3), grfullg(x, 2, 4)), c(0, 0), c(0, 0)),
cbind(c(0, 0, 0), c(0, 0, 0), c(0, 0, 0), 
 rbind(grfullg(x, 3, 5), grfullg(x, 3, 6), grfullg(x, 3, 7)))
)


print(resjacphi[1:n, (n+1):(n+m)] - checkB)	


cat("\n\n________________________________________\n\n")
cat("\n\npart C\n\n")	


gx <- c(g(x,1), g(x,2), g(x,3))

checkC <- 
- t(
cbind(
rbind(
grfullg(x, 1, 1) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 2) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 3) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 4) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 5) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 6) * GrAphiFB(-gx, lam)[1],
grfullg(x, 1, 7) * GrAphiFB(-gx, lam)[1]
),
rbind(
grfullg(x, 2, 1) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 2) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 3) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 4) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 5) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 6) * GrAphiFB(-gx, lam)[2:3],
grfullg(x, 2, 7) * GrAphiFB(-gx, lam)[2:3]
),
rbind(
grfullg(x, 3, 1) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 2) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 3) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 4) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 5) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 6) * GrAphiFB(-gx, lam)[4:5],
grfullg(x, 3, 7) * GrAphiFB(-gx, lam)[4:5]
)
)
)



print(resjacphi[(n+1):(n+m), 1:n] - checkC)


cat("\n\n________________________________________\n\n")

cat("\n\npart D\n\n")	


checkD <- diag(GrBphiFB(-gx, lam)) 

print(resjacphi[(n+1):(n+m), (n+1):(n+m)] - checkD)

Run the code above in your browser using DataLab