search

MOSEK / Tutorials

A collection of tutorials for the MOSEK package
107 stars 31 forks source link

Bug of Exponential cone? --- MORE --- src/prslv/checkpro.c:6177 #7

Closed hzzzzjzyq closed 5 months ago

hzzzzjzyq commented 5 months ago

Question: When we run the solution using mosek, two of the 10 calculation examples cannot be solved, and the error is --- MORE --- src/prslv/checkpro.c:6177. Can you help me solve it? What causes this? Is this a Mosek bug? Thank you for your help

''' import numpy as np import mosek import sys, time, math

def read_ins(filename): with open(filename, mode='r') as f: lines = [line.rstrip().split() for line in f.readlines() if len(line) > 1] lines.pop(0)

m, n = len(lines[1])//2, len(lines)

m,n=c_su,num_accept
a = [[float(lines[j][2*i]) for j in range(n)] for i in range(m)]
b = [[float(lines[j][2*i+1]) for j in range(n)] for i in range(m)]
return a, b

scale=100 c_su = 1000 num_accept = 10

r_list=[3.0,6.0,9.0,12.0,15.0]

r_list=[3.0]

file = open('zzzzz1224' + str(c_su) + '_print' + '.txt', 'a') for r in r_list: for k in range(10): tt=k+1 if tt==10: a, b = read_ins( '/home/zzjy/zzjy_file/pycharm_project/Kaggle-American-Express-Default-Prediction-1st-solution/instance/ins_large/data_1000010' + str(tt)) else: a, b = read_ins('/home/zzjy/zzjy_file/pycharm_project/Kaggle-American-Express-Default-Prediction-1st-solution/instance/ins_large/data_1000_010_0'+str(tt))

a, b = read_ins('/home/zzjy/zzjy_file/pycharm_project/Kaggle-American-Express-Default-Prediction-1st-solution/instance/ins_small/data_0050_005_01')

    a = list(map(list, zip(*a)))
    b = list(map(list, zip(*b)))

    # 将通过率和坏账率矩阵转换为numpy数组
    a_matrix = np.asarray(a)
    b_matrix = np.asarray(b)

    #本金
    f=1000000
    #利息收入率
    g=0.08

    ##############
    # 创建Mosek模型
    from mosek.fusion import *

    import signal
    import time

    def handler(signum, frame):
        raise Exception("Timeout")

    # 设置超时时间为10秒
    signal.signal(signal.SIGALRM, handler)
    signal.alarm(12)

    try:
        with Model() as m:
            # 定义变量
            m.setSolverParam("optimizerMaxTime", 10.0)
            m.setLogHandler(sys.stdout)  # 将日志输出到控制台

            x = [[m.variable('x_' + str(t) + '_' + str(c), Domain.binary()) for c in range(c_su)] for t in
                 range(num_accept)]

            obj = m.variable("obj", Domain.unbounded())  # Ensure z is greater than or equal to 0

            # 指数锥1
            z = m.variable("z", Domain.unbounded())  # Ensure z is greater than or equal to 0

            sum_term3 = [Expr.mul(-(g + 1) / (r * g)*scale, Expr.mul(b_matrix[t][c], x[t][c])) for t in range(num_accept) for
                         c
                         in range(c_su)]
            sum_expr3 = Expr.add(1*scale, Expr.add(sum_term3))  # 单独解第一项，用锥解的好一点点，没有小数后的误差。
            m.constraint('cone', Expr.vstack(sum_expr3, 1, z), Domain.inPExpCone())

            # 第二项
            objective_terms_2 = [
                Expr.mul(math.log(a_matrix[t][c]), x[t][c]) for t in range(num_accept) for c in range(c_su)
            ]

            # 指数锥2
            term_e = Expr.add(objective_terms_2)

            # 目标函数1
            m.objective(ObjectiveSense.Maximize, Expr.add(term_e, z))

            # 添加约束，确保每列的变量之和小于等于 1
            for c in range(c_su):
                sum_term6 = [x[t][c] for t in range(num_accept)]
                m.constraint(Expr.add(sum_term6), Domain.lessThan(1.0))

            # 添加约束，确保 sum_term2 之和小于等于 1
            sum_term4 = [Expr.mul((g + 1) / (r * g), Expr.mul(b_matrix[t][c], x[t][c])) for t in range(num_accept) for c
                         in range(c_su)]
            m.constraint(Expr.add(sum_term4), Domain.lessThan(1.0))

            sum_term5 = [x[t][c] for t in range(num_accept) for c in range(c_su)]
            m.constraint(Expr.add(sum_term5), Domain.equalsTo(r))
            # 求解

            start_time = time.time()
            m.solve()
            ourtime = time.time() - start_time
            print(ourtime)

            m.writeTask('path_to_output_file' + str(time.time()) + '.ptf')

            # 检查问题状态
            status = m.getProblemStatus()
            sol = []
            if status == mosek.fusion.ProblemStatus.PrimalInfeasible:
                print("原始问题不可行")
            else:
                for t in range(num_accept):
                    for c in range(c_su):
                        if x[t][c].level() > 0.5:
                            print(f'x_{t}_{c} =', x[t][c].level())
                            sol.append((t, c))
                obj_value = m.primalObjValue

                print("目标函数值:", obj_value)
        print('ok')

        product = 1
        for cusol in sol:
            product *= a_matrix[cusol[0]][cusol[1]]
        objvalue = f * g * (1 - (g + 1) / (r * g) * sum(b_matrix[cusol[0], cusol[1]] for cusol in sol)) * product

        print(objvalue)
        pass
    except Exception as e:
        print(f"Code execution took too long: {e}")
        objvalue=1
        ourtime=1
    finally:
        # 重置闹钟
        signal.alarm(0)

print('zzzzzc_su:'+str(c_su)+'_r:'+str(r)+'_num_accept:'+str(num_accept)+'_k:'+str(k)+'obj:'+str(objvalue)+'_time:'+str(ourtime)+'_scale:'+str(scale),file=file) print('---------------------------end----------------------------') file.close() ''' Problem Name: Objective sense: maximize Type: CONIC (conic optimization problem) Constraints: 1002 Affine conic cons. : 1 (3 rows) Disjunctive cons. : 0 Cones: 0 Scalar variables: 10003 Matrix variables: 0 Integer variables: 10000

Optimizer started. Mixed integer optimizer started. Threads used: 6 Presolve started. Presolve terminated. Time = 0.88, probing time = 0.81 Presolved problem: 9812 variables, 1002 constraints, 29428 non-zeros Presolved problem: 0 general integer, 9809 binary, 3 continuous Presolved problem: 1 cones Clique table size: 1038 BRANCHES RELAXS ACT_NDS DEPTH BEST_INT_OBJ BEST_RELAX_OBJ REL_GAP(%) TIME 0 1 1 0 NA -1.5246018912e+00 NA 0.9 --- MORE --- src/prslv/checkpro.c:6177 Code execution took too long: _PyEval_EvalFrameDefault returned a result with an error set exp_cone

  1. We transform problem (12) into an exponential cone problem. Where g, a, and b are the given data, g=0.08, a and b are numbers between 0 and 1 in 10 rows and 1000 columns.
  2. Data set: We generated 10 calculation examples, two of which could not be solved.
  3. We have tried a method: Since the term log(1-(g+1)/g bx) may be too small to cause numerical problems, we multiply the value in log by scale=100. But this doesn't completely solve the problem.
erling-d-andersen commented 5 months ago

That is a bug. Are you using the latest Mosek version?

søn. 19. maj 2024 04.15 skrev hzzzzjzyq @.***>:

Question: When we run the solution using mosek, two of the 10 calculation examples cannot be solved, and the error is --- MORE --- src/prslv/checkpro.c:6177. Can you help me solve it? What causes this? Is this a Mosek bug? Thank you for your help

''' import numpy as np import mosek import sys, time, math

def read_ins(filename): with open(filename, mode='r') as f: lines = [line.rstrip().split() for line in f.readlines() if len(line) > 1] lines.pop(0)

m, n = len(lines[1])//2, len(lines)

m,n=c_su,num_accept a = [[float(lines[j][2 i]) for j in range(n)] for i in range(m)] b = [[float(lines[j][2i+1]) for j in range(n)] for i in range(m)] return a, b

scale=100 c_su = 1000 num_accept = 10 r_list=[3.0,6.0,9.0,12.0,15.0]

r_list=[3.0]

file = open('zzzzz1224' + str(c_su) + '

print' + '.txt', 'a') for r in r_list: for k in range(10): tt=k+1 if tt==10: a, b = read_ins( '/home/zzjy/zzjy_file/pycharm_project/Kaggle-American-Express-Default-Prediction-1st-solution/instance/ins_large/data_1000_010'

  • str(tt)) else: a, b = read_ins('/home/zzjy/zzjy_file/pycharm_project/Kaggle-American-Express-Default-Prediction-1st-solution/instance/ins_large/data_1000_010_0'+str(tt))

    a, b =

    read_ins('/home/zzjy/zzjy_file/pycharm_project/Kaggle-American-Express-Default-Prediction-1st-solution/instance/ins_small/data_0050_005_01')

    a = list(map(list, zip(a))) b = list(map(list, zip(b)))

    将通过率和坏账率矩阵转换为numpy数组

    a_matrix = np.asarray(a) b_matrix = np.asarray(b)

    本金

    f=1000000

    利息收入率

    g=0.08

    ##############

    创建Mosek模型

    from mosek.fusion import *

    import signal import time

    def handler(signum, frame): raise Exception("Timeout")

    设置超时时间为10秒

    signal.signal(signal.SIGALRM, handler) signal.alarm(12)

    try: with Model() as m:

    定义变量

        m.setSolverParam("optimizerMaxTime", 10.0)
    m.setLogHandler(sys.stdout)  # 将日志输出到控制台

    x = [[m.variable('x_' + str(t) + '_' + str(c), Domain.binary()) for c in range(c_su)] for t in
         range(num_accept)]

    obj = m.variable("obj", Domain.unbounded())  # Ensure z is greater than or equal to 0

    # 指数锥1
    z = m.variable("z", Domain.unbounded())  # Ensure z is greater than or equal to 0

    sum_term3 = [Expr.mul(-(g + 1) / (r * g)*scale, Expr.mul(b_matrix[t][c], x[t][c])) for t in range(num_accept) for
                 c
                 in range(c_su)]
    sum_expr3 = Expr.add(1*scale, Expr.add(sum_term3))  # 单独解第一项，用锥解的好一点点，没有小数后的误差。
    m.constraint('cone', Expr.vstack(sum_expr3, 1, z), Domain.inPExpCone())

    # 第二项
    objective_terms_2 = [
        Expr.mul(math.log(a_matrix[t][c]), x[t][c]) for t in range(num_accept) for c in range(c_su)
    ]

    # 指数锥2
    term_e = Expr.add(objective_terms_2)

    # 目标函数1
    m.objective(ObjectiveSense.Maximize, Expr.add(term_e, z))

    # 添加约束，确保每列的变量之和小于等于 1
    for c in range(c_su):
        sum_term6 = [x[t][c] for t in range(num_accept)]
        m.constraint(Expr.add(sum_term6), Domain.lessThan(1.0))

    # 添加约束，确保 sum_term2 之和小于等于 1
    sum_term4 = [Expr.mul((g + 1) / (r * g), Expr.mul(b_matrix[t][c], x[t][c])) for t in range(num_accept) for c
                 in range(c_su)]
    m.constraint(Expr.add(sum_term4), Domain.lessThan(1.0))

    sum_term5 = [x[t][c] for t in range(num_accept) for c in range(c_su)]
    m.constraint(Expr.add(sum_term5), Domain.equalsTo(r))
    # 求解

    start_time = time.time()
    m.solve()
    ourtime = time.time() - start_time
    print(ourtime)

    m.writeTask('path_to_output_file' + str(time.time()) + '.ptf')

    # 检查问题状态
    status = m.getProblemStatus()
    sol = []
    if status == mosek.fusion.ProblemStatus.PrimalInfeasible:
        print("原始问题不可行")
    else:
        for t in range(num_accept):
            for c in range(c_su):
                if x[t][c].level() > 0.5:
                    print(f'x_{t}_{c} =', x[t][c].level())
                    sol.append((t, c))
        obj_value = m.primalObjValue

        print("目标函数值:", obj_value)
print('ok')

product = 1
for cusol in sol:
    product *= a_matrix[cusol[0]][cusol[1]]
objvalue = f * g * (1 - (g + 1) / (r * g) * sum(b_matrix[cusol[0], cusol[1]] for cusol in sol)) * product

print(objvalue)
pass

    except Exception as e: print(f"Code execution took too long: {e}") objvalue=1 ourtime=1 finally:

    重置闹钟

    signal.alarm(0)

print('zzzzzc_su:'+str(c_su)+'_r:'+str(r)+'_num_accept:'+str(num_accept)+'_k:'+str(k)+'obj:'+str(objvalue)+'_time:'+str(ourtime)+'_scale:'+str(scale),file=file) print('---------------------------end----------------------------') file.close() ''' Problem Name: Objective sense: maximize Type: CONIC (conic optimization problem) Constraints: 1002 Affine conic cons. : 1 (3 rows) Disjunctive cons. : 0 Cones: 0 Scalar variables: 10003 Matrix variables: 0 Integer variables: 10000

Optimizer started. Mixed integer optimizer started. Threads used: 6 Presolve started. Presolve terminated. Time = 0.88, probing time = 0.81 Presolved problem: 9812 variables, 1002 constraints, 29428 non-zeros Presolved problem: 0 general integer, 9809 binary, 3 continuous Presolved problem: 1 cones Clique table size: 1038 BRANCHES RELAXS ACT_NDS DEPTH BEST_INT_OBJ BEST_RELAX_OBJ REL_GAP(%) TIME 0 1 1 0 NA -1.5246018912e+00 NA 0.9 --- MORE --- src/prslv/checkpro.c:6177 Code execution took too long: _PyEval_EvalFrameDefault returned a result with an error set exp_cone.png (view on web) https://github.com/MOSEK/Tutorials/assets/45667990/ee0489e0-a83e-4f42-bd75-3d4dcb654f2a

  1. We transform problem (12) into an exponential cone problem. Where g, a, and b are the given data, g=0.08, a and b are numbers between 0 and 1 in 10 rows and 1000 columns.
  2. Data set: We generated 10 calculation examples, two of which could not be solved.
  3. We have tried a method: Since the term log(1-(g+1)/g bx) may be too small to cause numerical problems, we multiply the value in log by scale=100. But this doesn't completely solve the problem.

— Reply to this email directly, view it on GitHub https://github.com/MOSEK/Tutorials/issues/7, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACOAE4PQW7L6TNFSY5ORD43ZDADM5AVCNFSM6AAAAABH56KGFSVHI2DSMVQWIX3LMV43ASLTON2WKOZSGMYDIMZZGAYTOOI . You are receiving this because you are subscribed to this thread.Message ID: @.***>

erling-d-andersen commented 5 months ago

To me it seems this problem should not occur if you use MOSEK version 10.1.36.

If you think the problem is still there, then please dump the problem to disk using the instructions

https://docs.mosek.com/latest/faq/faq.html

email it to support@mosek.com.

erling-d-andersen commented 5 months ago

I close the issue since I have not heard anything.