From e7c649964dae3cbc2e41bfe5d212d093ef8e006c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= Date: Sat, 17 Jan 2026 16:13:17 +0100 Subject: [PATCH 1/2] Soft constraints for AVIs --- interfaces/daqp-julia/src/avi.jl | 115 +++++++++++++++---------------- 1 file changed, 54 insertions(+), 61 deletions(-) diff --git a/interfaces/daqp-julia/src/avi.jl b/interfaces/daqp-julia/src/avi.jl index 2a19ca9..b83c59a 100644 --- a/interfaces/daqp-julia/src/avi.jl +++ b/interfaces/daqp-julia/src/avi.jl @@ -10,16 +10,20 @@ end struct AVIWorkspace H::Matrix{Float64} f::Vector{Float64} - A::Matrix{Float64} + At::Matrix{Float64} bu::Vector{Float64} bl::Vector{Float64} + sense::Vector{Cint} + norm_factors::Vector{Cdouble} - H_factor::LU{Float64} + rho_soft::Float64 + + H_factor::LU{Float64, Matrix{Float64}, Vector{Int64}} x_unc::Vector{Float64} H1pI::Matrix{Float64} H2mI::Matrix{Float64} - H2pI::LU{Float64} # LU + H2pI::LU{Float64, Matrix{Float64}, Vector{Int64}} H2::Matrix{Float64} H2xpf::Vector{Float64} x::Vector{Float64} @@ -30,13 +34,13 @@ struct AVIWorkspace settings::AVISettings end -AVIWorkspace() = AVIWorkspace(zeros(0,0),zeros(0),zeros(0,0),zeros(0),zeros(0), - lu(zeros(0,0)),zeros(0), +AVIWorkspace() = AVIWorkspace(zeros(0,0),zeros(0),zeros(0,0),zeros(0),zeros(0),zeros(Cint,0),zeros(0), + 0,lu(zeros(0,0)),zeros(0), zeros(0,0),zeros(0,0),lu(zeros(0,0)),zeros(0,0), zeros(0), zeros(0),zeros(0), Model(), zeros(0), AVISettings()) -function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), daqp_workspace=nothing, settings=AVISettings()) +function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), rho_soft = 1e6, daqp_workspace=nothing, settings=AVISettings()) n,m = length(f),length(bu) bl = isempty(bl) ? fill(-1e30,m) : bl sense = isempty(sense) ? zeros(Cint,m) : sense @@ -51,8 +55,8 @@ function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), daqp_workspace=nothing, sett x = isempty(x0) ? fill(NaN,n) : copy(x0) y,H2xpf = zeros(n),zeros(n) - kkt_buffer = Float64[] - sizehint!(kkt_buffer,4*(n+1)*(n+1)) + kkt_buffer = zeros(4*(n+1)*(n+1)) + #sizehint!(kkt_buffer,4*(n+1)*(n+1)) daqp_workspace = isnothing(daqp_workspace) ? Model() : daqp_workspace; DAQPBase.settings(daqp_workspace, Dict(:iter_limit => settings.inner_iter_limit, @@ -60,17 +64,43 @@ function setup_avi(H,f,A,bu,bl,sense;x0 = zeros(0), daqp_workspace=nothing, sett :dual_tol => settings.dual_tol)) exitflag,tsetup = DAQPBase.setup(daqp_workspace, H1pI, H2xpf, A, bu, bl, sense) + workspace = unsafe_load(daqp_workspace.work); + norm_factors = unsafe_wrap(Vector{Cdouble}, workspace.scaling, m; own = false) + + if m > size(A,1) - A = [I(n)[1:m-size(A,1),:];A] # To handle simple bounds + #A = [I(n)[1:m-size(A,1),:];A] # To handle simple bounds + At = [I(n)[:,1:m-size(A,1)] A'] # To handle simple bounds + else + At = A'[:,:] end - return exitflag, AVIWorkspace(H,f,A,bu,bl, - H_factor,zeros(n), + return exitflag, AVIWorkspace(H,f,At,bu,bl,sense,norm_factors, + rho_soft,H_factor,zeros(n), H1pI,H2mI,H2pI,H2,H2xpf, x,y,daqp_workspace,kkt_buffer, settings) end +function _is_optimal(λ,xt,λfull,ws,ASu,ASl) + nu = length(ASu) + ϵp, ϵd = ws.settings.primal_tol, ws.settings.dual_tol + @inbounds for i in 1:nu + λ[i] < -ϵd && return false + end + @inbounds for i in nu+1:length(λ) + λ[i] > ϵd && return false + end + @inbounds for i in 1:length(ws.bu) + if λfull[i] == 0.0 # Only test inactive constraints + Axi = dot(view(ws.At,:,i),xt) + Axi-ws.bu[i] > ϵp && return false + Axi-ws.bl[i] < -ϵp && return false + end + end + return true +end + function solve(ws::AVIWorkspace) exitflag = 0 λstar, ASstar = zeros(0), Int[] @@ -91,17 +121,12 @@ function solve(ws::AVIWorkspace) if info.iterations == 1 nkkt += 1 ASu,ASl = _get_AS(info.λ,ϵd) - #xt,λ,AS = _solve_kkt_schur(ws,ASu,ASl) xt,λ,AS = _solve_kkt(ws,ASu,ASl) - nu = length(ASu) - if all(λ[i] > -ϵd for i in 1:nu) && all(λ[i] < ϵd for i in nu+1:length(λ)) - Ax = ws.A*xt - if all(Ax-ws.bu .<= ϵp) && all(Ax-ws.bl .>=-ϵp) - ws.x .= xt - λstar,ASstar = copy(λ),AS - exitflag, outer_iter = 1,k - break - end + if _is_optimal(λ,xt,info.λ,ws,ASu,ASl) + ws.x .= xt + λstar,ASstar = copy(λ),AS + exitflag, outer_iter = 1,k + break end axpby!(β,xt,1.0-β,ws.x) else @@ -127,53 +152,21 @@ function _get_AS(λ,dual_tol) return ASu,ASl end -# TODO kkt_schur is faster in isolation, but leads to slower total execution -#function _solve_kkt_schur(ws::AVIWorkspace,ASu,ASl) -# # Prepare some helpers -# AS = [ASu;ASl] -# n,nAS = length(ws.f),length(AS) -# -# # Juggle some buffers -# resize!(ws.kkt_buffer,nAS*(nAS+n+1)+n) # Is this necessary? -# invH_At = reshape(view(ws.kkt_buffer,1:nAS*n),(n,nAS)) -# offset = nAS*n -# S = reshape(view(ws.kkt_buffer,offset+1:offset+nAS^2),(nAS,nAS)) -# offset +=nAS^2 -# λ = view(ws.kkt_buffer,offset+1:offset+nAS) -# offset +=nAS -# x = view(ws.kkt_buffer,offset+1:offset+n) -# -# AAS = ws.A[AS,:] -# -# # Schur complement -# @inbounds for (i,id) in enumerate(AS) -# @views invH_At[:,i] .= ws.A[id,:] -# end -# ldiv!(ws.H_factor,invH_At) -# mul!(S, AAS, invH_At) -# S_fact = lu!(S) -# -# # Solve for λ -# @views λ .= AAS * ws.x_unc -# @views λ[1:length(ASu)] .-= ws.bu[ASu] -# @views λ[length(ASu)+1:end] .-= ws.bl[ASl] -# ldiv!(S_fact,λ) -# -# # Solve for x -# x .= ws.x_unc .- invH_At * λ -# -# return x,λ,AS -#end function _solve_kkt(ws,ASu,ASl) - AS = [ASu;ASl] + AS = Int[ASu;ASl] n = length(ws.f) nkkt = n+length(AS) resize!(ws.kkt_buffer,nkkt^2) K = reshape(view(ws.kkt_buffer,1:nkkt^2),(nkkt,nkkt)) @views K[1:n,1:n] = ws.H - @views K[n+1:end,1:n] = ws.A[AS,:] - @views K[1:n,n+1:end] = ws.A[AS,:]' + @views K[n+1:end,1:n] = ws.At[:,AS]' + @views K[1:n,n+1:end] = ws.At[:,AS] K[n+1:end,n+1:end] .= 0 + for (i,id) in enumerate(AS) + if ws.sense[id] & DAQPBase.SOFT != 0 + K[n+i,n+i] = -1/((ws.norm_factors[id]^2)*ws.rho_soft) + end + end z = K\[-ws.f;ws.bu[ASu];ws.bl[ASl]] x = @view z[1:n] λ = @view z[n+1:end] From 93574ddec4debda13802d88cbc329cd8dac8312a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Daniel=20Arnstr=C3=B6m?= Date: Sat, 17 Jan 2026 16:13:58 +0100 Subject: [PATCH 2/2] Bump version --- interfaces/daqp-julia/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interfaces/daqp-julia/Project.toml b/interfaces/daqp-julia/Project.toml index 80a1148..ac853be 100644 --- a/interfaces/daqp-julia/Project.toml +++ b/interfaces/daqp-julia/Project.toml @@ -1,7 +1,7 @@ name = "DAQPBase" uuid = "833476c3-a8c0-4073-9b64-6473509843fe" authors = ["Daniel Arnström "] -version = "0.3.5" +version = "0.3.6" [deps] DAQP_jll = "5c51c916-43bf-57fe-9d62-13064ebbf40d"